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1.  INTRODUCTION  AND  SUMMARY 


The  objective  of  the  Multi  target  Tracking  Studies  (MTS)  project  is  to 
develop  and  evaluate  a  parameter  modeling  concept  for  tracking  multiple 
targets.  This  report  simmarizes  the  acconplishments  on  phase  II  of  the 
project,  during  the  period  of  June  1980  to  June  1981.  The  results  obtained 
during  phase  I  of  MTS  were  reported  in  [1]. 

The  MTS  concept  Is  based  on  modeling  the  observed  data  as  a  multichannel 
autoregressive  moving-average  (ARMA)  process.  The  parameters  of  the  model 
provide  a  compact  representation  of  target  parameters  such  as  spectrum  and 
time-difference-of-arrival  (TDOA).  These  parameters  can  be  used  as  inputs  to 
a  tracking  algorithm,  a  target  classification  program,  etc.  The  unique 
feature  of  the  MTS  technique  is  the  simultaneous  estimation  of  multiple  target 
parameters,  'lost  conventional  processing  techniques  estimate  parameters  of 
one  target  at  a  time,  considering  the  other  targets  as  sources  of 
interference. 

Parametric  modeling  of  ARMA  signals  at  low  signal-to-noise  ratios 
presents  a  difficult  and  challenging  problem,  even  in  the  single  channel 
case.  It  was  clear  from  the  start  that  single-input  single-output  ARMA  models 
have  to  be  studied  before  tackling  the  multi -Input  multi -output  models  that 
arise  in  the  multi  target  case.  Accordingly,  most  of  the  effort  in  phase  I  and 
part  of  the  effort  in  phase  II  was  devoted  to  the  single  channel  case.  The 
emphasis  was  on  developing  a  recursive  parameter  estimation  algorithm  capable 
of  working  at  the  low-signal -to-nolse  ratios  that  are  typically  encountered  In 
the  undersea  surveillance  scenario  (e.g.  -10  to  -15  dB).  In  section  2  we 
present  the  single  channel  algorithm  which  is  a  more  robust  version  of  the  one 
presented  in  [1].  This  algorithm  is  used  as  an  adaptive  Infinite  Impulse 
Response  (IIR)  filter  to  perform  a  minber  of  functions  such  as:  adaptive  line 
enhancement,  adaptive  noise  cancelling,  adaptive  TDOA  estimation,  adaptive 
deconvolution  and  high  resolution  spectral  estimation.  Adaptive  filters  of 
this  kind  are  useful  for  line  tracking,  own  ship  noise  reduction,  inter-array 
processing  and  other  surveillance  applications.  The  IIR  adaptive  filter  has 
Improved  performance  compared  to  conventional  Finite  Impulse 
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filter,  especially  under  low  slgnal-to-noise  conditions,  as  demonstrated  by 
the  results  In  [2]-[4]. 

In  section  3  we  describe  the  multichannel  ARMA  modeling  algorithm 
developed  for  the  multitarget  case.  The  algorithm  was  tested  under  high 
signal-to-noise  conditions  and  was  shown  to  be  able  to  estimate  correctly  the 
TDOA's  and  spectral  parameters  of  two  targets.  However,  a  number  of 
theoretical  as  well  as  practical  issues  remain  to  be  resolved  before  the 
performance  of  the  algorithm  can  be  adequately  assessed.  This  will  be  done  as 
part  of  phase  III  of  the  MTS  project  in  which  a  more  robust  version  of  the 
multi  target  processor  will  be  developed.  The  main  difficulties  encountered  in 
the  ARMA  approach  are  related  to  the  question  of  structural  constraints.  The 
parametric  tracking  model  (comprised  of  a  spectral  model  and  a  propagation 
model)  considered  in  our  work  leads  to  a  certain  linear  transfer  function 
matrix  describing  the  signals  received  by  the  sensors.  The  transfer  matrix  is 
called  the  Right  Matrix  Fraction  Description  (RMFD)  of  the  system.  The  RMFD 
associated  with  the  tracking  model  has  a  very  special  structure  as  will  be 
discussed  in  section  3.  Estimating  the  parameters  of  this  matrix  involves  the 
solution  of  a  highly  nonlinear  optimization  problem.  The  4RMA  model  provides 
an  alternative  parameterization  (a  Left  Matrix  Fraction  Description)  which  is 
much  easier  to  compute.  Unfortunately,  the  special  structure  of  the  RMFD  is 
not  shared  by  the  ARMA  model.  The  algorithm  we  currently  use  estimates  the 
parameters  of  a  general  ARMA  model  wfth  no  structural  constraints.  Because  of 
this,  the  tracking  model  is  highly  overparametrized,  leading  to  uniqueness 
problems  and  to  degraded  performance  at  low  SNR.  To  make  the  MTS  processor 
more  robust  it  seems  that  we  will  have  to  replace  the  recursive  parameter 
estimation  algorithm  originally  used  with  a  different  type  of  algorithm 
capable  of  enforcing  the  special  model  structure.  This  will  probably  lead  to 
a  more  complex  algorithm  structure. 

The  research  performed  so  far  on  the  MTS  project  resulted  in  a  number  of 
publications,  as  listed  in  section  5.  In  this  report  we  present  a  brief 
sunmary  of  this  research,  deferring  many  details  to  these  publications.  Two 
of  the  key  papers,  one  for  the  single  channel  case  and  one  for  the  multi  target 
case,  are  included  as  appendices  A  and  B.  Further  details  on  the  MTS 
algorithms  are  presented  in  appendices  C  and  D. 
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2.  SINGLE  TARGET  ALGORITHMS 
2.1  THE  RML  ALGORITHM 

The  MTS  processor  is  based  on  fitting  an  ARMA  model  to  an  observed  time- 
series.  Several  recursive  parameter  estimation  algorithms  were  studied  in 
phase  I  for  performing  this  model  fitting.  The  recursive  maximum  likelihood 
algorithm  (RML)  was  chosen  as  the  one  most  suitable  for  our  problem.  The  RML 
algorithm  estimates  the  parameters  of  an  ARMA  model  of  the  type 

NA  MB  NC 

h  ■  -£  a<  Vi  *  H  bi  Vi  *  £  ci  Vi  *  vt  m 

i=l  i=l  i=l 


where  {y^,  u^}  are  the  data  sequence  and  vt  is  an  unobservable  white  noise 
sequence.  The  RML  algorithm  is  specified  by  the  following  set  of  equations 
(see  Appendix  A) 

0*  [a1,...,aMA,  b1,...,bNB,  c1 . c^c]T  *  parameter  vector  (2a) 

*t  *  t'yt-l’"-,“yt-NA’  ut-l*  * '  *  ,ut-NB’  et-r- **  ^t-NC-1  s  vector  (2bl 

\  *  ^~yt-l  ’ ’ ‘  ’  ’~yt-NA,  “t-l*"*’  Ut-NB *  et-l . et-NC-*  "  data  vector  2  ' 

pt  *  *-Pt-l  "  Pt-lli;t  ‘H;  Pt-l^xt  +  ^t  Pt-l(t't'^xt 

9t  =  9t-l  +  pt  h  (yt  "  *t  Vl’  (4) 


et  =  yt  •  h  9t 


yt  =  (1/Dt(z))  yt,  Zt  =  (1/Dt(z))  yt,  et  =  (l/Ot(z) )  et 
The  pre-filter  Dt(z)  is  usually  taken  to  be 
Dt(z)  ■  1+c^t)  z-1+...+  cNC(t)  z'NC 


5 


The  expontential  forgetting  factor  is  usually  chosen  as, 
Xt  3  X  \-l  + 

and  the  initial  conditions  are 


Several  features  were  added  to  the  algorithm  to  guarantee  its  convergence 
and  make  it  more  robust: 

(i)  Stability  monitoring 

The  stability  of  the  pre-filter  Dt(z)  needs  to  be  tested.  Whenever  Dt(z) 
becomes  unstable  the  parameter  estimates  9t  need  to  be  projected  back  into  a 
region  of  stability  [5]. This  can  be  done  by  setting  the  parameters  back  to 
their  value  before  the  last  update:  e*  =  e«.  ,  Another  procedure  is  to  go 
back  only  part  of  the  way,  e.g. , 


4t  *  at-i  -sVt 


where 


vt  *  \-i  *t  (yt  -  *t  Vi1 


(8a) 

(8b) 


The  coefficient  s  is  set  initially  to  max  {1/2,  10'6/|Vtl),  where  | V^l 
denotes  the  sum  of  the  absolute  values  of  the  entries  of  the  vector  V^..  If 
the  new  parameter  estimates  9*  are  still  unstable  we  set 

£  V 

s  =  max  {s/2,  10  /|VJ)  and  repeat  the  process.  The  stability  test 

t 

algorithm  is  given  by  the  following  program: 
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FUNCTION  NSTAB(T.NT) 


C  PERFORMS  A  TEST  WHETHER  THE  POLYNOMIAL  T  HAS  ALL  ROOTS 
C  INSIDE  THE  UNIT  CIRCLE  OR  NOT 
C  T  *  1+T(1)*Q**(-1)+. . .+T(NT)*Q**{-NT) 

C  WHERE  Q** ( - 1 )  IS  THE  BACKWARD  SHIFT  OPERATOR 
C 

C  T  IS  A  VECTOR  CONTAINING  THE  COEFFICIENTS  OF  THE 
C  POLYNOMIAL 

C  NT  IS  THE  ORDER  OF  THE  POLYNOMIAL  /MAX=30/ 

C 

C  THE  FUNCTION  IS  RETURNED  0  IF  ALL  ROOTS  ARE  INSIDE  THE 
C  UNIT  CIRCLE,  I.E.  FOR  A  STABLE  POLYNOMIAL;  OTHERWISE  AS  1 
DIMENSION  T(l) 

DIMENSION  TT(30) ,  TTT(30) 

C 

NSTAB=0 

r 

DO  10  1=1,  NT 
10  TT ( I ) =T ( I ) 

C 

00  50  I=NT,1 , -1 
AK=TT ( I ) 

IF(ABS(AK).GE.l.)  GO  TO  60 
Tl*l. -AK*AK 
11=1-1 

IF(I1)  50,50,20 
20  DO  30  J=1,I1 

30  TTT(J)=(TT(J)-AK*TT(I-J) )/Tl 
DO  40  J-1,11 
40  TT(J )  =  TTT(J) 

50  CONTINUE 

C 

GO  TO  70 
60  NSTAB*1 
70  RETURN 
END 
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( i  1 )  Modified  Pre-filtering 


The  choice  of  the  pre-filter  Dt(z)  has  an  effect  both  on  the  asymptotic 
properties  of  the  RML  algorithm  and  on  its  transient  behavior.  The  choice 
indicated  in  equation  (7)  guarantees  asymptotic  convergence  and  asymptatic 
efficiency.  However,  other  choices  seem  to  have  faster  convergence  rate, 
which  is  an  important  factor  in  adaptive  processing  applications. 

Me  experimented  with  a  modified  pre-filter  of  the  form: 

Dt(z)  =  1  +  kc1(t)z'1+...+lcNC  cNC(t)z‘NC  (9) 

where  0<k<l.  For  different  choices  of  the  parameter  k  we  get  algorithms  with 
different  types  of  behavior.  In  particular,  k*l  gives  the  RM.  algorithm  while 
k=G  gives  the  extended  least-squares  algorithm.  When  using  the  algorithm  on 
narrowband  signals  in  noise  we  found  that  using  intermediate  values  of  k  often 
improved  the  performance  of  the  algorithm.  A  more  detailed  discussion  of  the 
modified  pre-filter  Is  given  in  [6]. 

(iii)  The  Square-Root  Form  of  the  RML 

The  basic  form  of  the  Rft  Involves  the  updating  of  the  covariance  matrix 
Pt  using  the  difference  Eq.(4).  This  type  of  equation  does  not  guarantee  the 
positive  definiteness  of  Pt  and  can  lead  to  numerical  problems.  It  is  well 
known  that  a  much  better  behaved  solution  to  least-squares  estimation  problems 
can  be  obtained  by  square-root  techniques  which  update  p^rather  than 
Pt  (where  P^P^2  3  P^)  [7], [8].  We  used  the  following  algorithm  to 
Implement  a  square-root  version  of  the  RW. :  Let 

pt-l  *  U  0  *  old  covariance  matrix 

A  A  AT 

Pt  5  U  0  U  =  new  covariance  matrix 


where  O’,  u  are  upper  triangular  with  unity  on  the  diagonal,  and  d,  D  are 
diagonal  matrices.  Then  the  following  are  the  update  equations  for  the  ll-D 
factors: 


•  Initialization: 


O'  =  I,  D  =  diag  {a} 

•  First  step: 

f  =  UT  4»t.  where  fT  =  ffl] 

v  =  Of,  i.e.,  v.  =  d.  f.  for  i=l,...,  n 
=  d1/a^  where 
=  Cvlt  °,  ....  0] 

•  Ma i n  Loop 


for  j  =  2,  ....  n  do 


ao  *  *j-i  *  'j  fj 

dj =  dj  aj-i/uj x  t> 

“j  *"j  *aj  Vi'  '"’ere  aj  =  -fj/aj-i 
Ej  ■  EJ-1  *  vj  “J 

A  A 

We  use  the  notation  u^,  u^  for  the  columns  of  U,  U,  i.e., 
0”  »  [u^  ,  . ..,  un],  U  =  [Up  . ...  un] 

•  Conpute  New  Gain 


=  ^n^an 
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Update  Parameter  Vector 


et  =  9t-i  +  Kt  lyt  '  *t  9t-i* 

•  Compute  Residual 

et  *  ■  »t  9t 

An  implementation  of  these  update  equations  can  be  found  in  [8,  pp.  100- 
101].  A  slight  modification  was  needed  to  account  for  an  extra  division 
by  which  appears  in  (4),  but  not  in  the  standard  Kalman  update  of  the 
covariance  matrix  P^.. 

The  square-root  RM.  has  shown  somewhat  improved  performance,  especially 
at  low  SNR,  and  we  are  now  using  only  this  version  of  the  algorithm.  Since 
our  simulations  were  performed  using  high  precision  floating-point 
corrputuaions,  the  difference  between  the  regular  and  square-root  normalized 
RM.  is  not  large.  A  much  more  pronounced  differnce  is  expected  to  occur  when 
using  fixed  point  finite-word  length  computations. 

The  RM.  algorithm  with  the  various  featrues  described  above  is  quite 
robust  and  has  been  used  successfully  on  narrowband  data  at  SNR’s  as  low  as 
-15  dB.  The  algorithm  will  probably  work  at  lower  SNR  given  sufficient 
data.  In  our  simulations  we  used  not  more  that  4096  data  points  per  run. 

2.2  APPLICATIONS  OF  THE  RML  FOR  ADAPTIVE  PROCESSING 

The  RML  algorithm  can  be  used  for  a  variety  of  adaptive  processing 
applications.  In  appendix  A  we  present  a  summary  of  our  work  on  using  the  RML 
for  adaptive  line  enhancement,  adaptive  noise  cancelling,  adaptive  time  delay 
estimation,  adaptive  deconvolution  and  high  resolution  spectral  estimation. 
Here  we  provide  brief  descriptions  of  these  applications  and  give  references 
to  project  publications  containing  more  details. 

(i)  Adaptive  Line  Enhancement  [2] 
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The  RML  algorithm  can  be  used  as  an  adaptive  prediction  filter  for  ARMA 
processes.  The  predicted  value  y^^  is  given  by 
NA  NC 

i=l  i=l 

where 

T  : 

£t-i  "  yt-l  ■  Vrt-i-1  * 

The  structure  of  the  IIR  filter  is  depicted  in  Fig.  1.  The  coefficients 
of  this  filter  are  adjusted  at  each  time  step  according  to  the  update  equation 
(2)  -  (7)  of  the  RM.  algorithm.  When  the  input  is  a  narrowband  signal 
corrupted  by  white  noise  the  adaptive  predictor  will  form  a  bandpass  filter 
around  the  spectral  lines  and  reject  much  of  the  noise.  In  other  words,  the 
adaptive  predictor  is  precisely  the  so-called  adaptive  line  enhancer  (ALE) 

[9]. 


The  ALE's  currently  used  by  the  Navy  and  other  users  of  adaptive  signal 
processing  are  all  based  on  FIR  filtering.  As  discussed  in  [2],  the  IIR 
filter  has  a  potential  for  performance  improvement  over  FIR  filters, 
especially  at  low  SNR's.  The  SNR  improvement  of  the  FIR-ALE  is  proportional 
to  the  filter  order.  Thus,  to  achieve  significant  noise  rejection,  filters  of 
very  high  order  have  to  be  used,  leading  to  fairly  complex  hardware.  Filters 
with  4096  weights  and  more  are  not  uncommon.  The  IIR-ALE  can  achieve  better 
performance  with  a  very  low  order  filter,  the  order  being  just  twice  the 
nunber  of  spectral  lines  of  interest.  This  can  lead  to  significant 
computational  saving  and  simpler  hardware.  We  have  also  developed  a  lattice 
version  of  the  IIR-ALE  which  is  particularly  attractive  for  special  purpose 
hardware  implementation  [10]. 

( 1 i )  Adaptive  Noise  Cancelling  [3] 

The  noise  cancelling  technique  makes  use  of  an  auxiliary  or  reference 
Input  which  provides  some  information  about  the  noise  component  that  is 
corrupting  the  signal.  From  this  "side  information"  an  estimate  of  the  noise 
is  obtained  and  then  subtracted  from  the  primary  input  to  "cancel  out"  the 
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additive  noise  component.  The  reference  input  may  be  an  extra  beam  steered  in 
the  direction  of  known  noise  sources  (e.g.  own  ship)  or  extra  sensors  in  a 
towed  array  giving  additional  information  about  flow  noise.  These  noise 
subtraction  techniques  are  very  sensitive  to  modeling  errors  since  they  are 


'NPUT  £SR0* 


Figure  1.  The  IIR  Adaptive  Line  Enhancer. 


essentially  equivalent  to  "bridge  balancing".  Therefore,  noise  cancelling  is 
typically  performed  by  an  adaptive  processor  which  properly  adjusts  the  filter 
parameters. 


All  of  the  adaptive  noise  cancellers  (ANC)  currently  in  use  seem  to  be  of 
the  FIR  type.  Using  the  RMl  algorithm  we  were  able  to  develop  an  IIR-ANC 
whose  structure  is  depicted  in  Fig.  2.  In  [3]  it  was  shown  that  the  general 
ANC  problem  can  be  solved  by  fitting  the  following  model  to  the  observed  data 
at  the  primary  (yt)  and  reference  inputs  (ut): 


*t 


F(z)  r 


(11) 


where  et  i$  an  unmeasurable  white  noise  process.  A  modified  version  of  the 
RML  algorithm  [11]  recursively  estimates  the  parameters  of  this  model,  which 
are  then  used  to  adjust  the  filter  coefficients. 


Note  that  the  IIR-ANC  involves  only  the  parameters  of  B(z)  and  F ( z) -  It 
turns  out  that  the  remaining  parameters  C(z),  D(z)  can  be  used  to  form  an 
adaptive  line  enhancer  for  the  signal  at  the  ANC  ouput.  In  other  words  the 
RM.  algorithm  can  adjust  simultaneously  an  IIR-ANC  and  an  IIR-ALE,  as  depicted 
in  Fig.  3.  The  combined  ALE-ANC  was  tested  on  synthetic  data  with  promising 
results,  as  reported  in  [3], 
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combined  ALE-ANC. 


(iii)  High  Resolution  Spectral  Estimation  [4] 

Spectral  estimation  techniques  play  an  important  role  in  dectection, 
localization  and  classification  of  targets.  High  resolution  AR  and  ARMA 
methods  were  used  in  many  applications  such  as  seismic  data  processing,  radar 
and  sonar  processing  and  general  time-series  analysis  problems.  Several  good 
techniques  are  currently  available  for  AR  modeling  of  time  series.  The  ARMA 
modeling  problem  is  much  more  difficult,  and  leads  to  optimization  of  a  highly 
nonlinear  likelihood  function.  Because  of  the  computational  complexity  of  the 
exact  maximum  likelihood  technique,  ARMA  modeling  is  often  avoided  in  practice 
in  spite  of  its  potential  for  Improved  resolution  and  better  spectral 
matching.  Various  suboptimal  batch  processing  techniques  were  proposed  in  the 
past  for  ARMA  spectral  estimation.  The  RML  algorithm  provides  a  suboptimal 
recursive  solution  technique  which  is  well  suited  for  adaptive  processing 
applications.  We  tested  the  RML  spectral  estimator  on  several  types  of 
signals  under  various  SNR  conditions.  Results  have  been  good  for  sufficiently 
long  data  records,  as  reported  in  [4].  On  short  data  records  recursive 
techniques  are  bound  to  be  Inferior  to  batch  techniques,  because  of  the 
transient  behavior  of  such  algorithms.  The  RML  can  also  be  used  in  a  batch 
mode,  in  which  case  it  is  passed  several  times  over  a  given  data  record.  When 
used  in  this  mode  the  RML  spectral  estimates  can  match  those  of  competitive 
ARMA  techniques. 

Signal  processing  problems  often  involve  sinusoidal  signals  corrupted  by 
additive  noise.  When  it  Is  known  a  priori  that  the  signals  of  interest  are 
sinusoids,  this  Information  can  be  used  to  enhance  the  performance  of  spectral 
estimators.  When  AR  techniques  are  used, the  presence  of  sinusoids  means  that 
the  roots  of  the  predictor  A(z)  lie  on  the  unit  circle.  Another 
interpretation  of  this  property  is  that  A(z)  Is  a  linear  phase  filter. 

Forcing  the  AR  prediction  filter  to  have  linear  phase  characteristics  leads  to 
improved  spectral  estimation  performance  as  was  shown  by  Marple  [12].  In  [13] 
we  developed  several  recursive  least-squares  algorithms  for  AR  modeling  which 
can  be  used  for  spectral  estimation  and  other  applications. 


(iv)  Adaptive  Deconvolution 


The  need  to  extract  a  signal  from  a  filtered  version  of  that  signal 
arises  in  many  situations  such  as  channel  equalization,  seismic  data 
processing,  speech  analysis,  and  removal  of  shallow  water  reflections  in  an 
active  sonar  system.  The  RML  algorithm  can  be  directly  applied  to  these 
problems  provided  that  the  signal  (before  it  was  filtered)  is  sufficiently 
white.  In  this  case  the  RML  algorithm  is  used  to  implement  an  adaptive 
',whitener,,.  The  prediction  error  sequence  will  then  provide  an  estimate  of 
the  original  signal.  Note  that  this  technique  can  handle  signals  that  have 
passed  a  pole- zero  filter,  whereas  most  of  the  current  deconvolution 
techniques  are  limited  to  all-pole  filters.  Some  examples  of  adaptive 
deconvolution  using  the  RML  algorithm  can  be  found  in  appendix  A.  In  speech 
analysis  and  other  applications  it  is  known  a  priori  that  the  signal  is  a 
sequence  of  impulses  of  random  anplitudes  occuring  at  random  points  in  time. 
This  knowledge  can  be  used  to  improve  the  performance  of  the  deconvolution 
technique  by  performing  a  nonlinear  opertion  in  the  feedback  loop  within  the 
RML  algorithm;  see  [14]  for  details. 

(v)  Adaptive  Time  Delay  Estimation 

Delay  estimation  is  an  important  technique  for  target  localization  by 
multiple  sensors  or  multiple  arrays.  Many  different  techniques  for  delay 
estimation  have  been  proposed,  as  can  be  seen  In  [15].  These  techniques 
involve  pre-filtering  followed  by  cross  correlation  and  are  usually  performed 
in  the  frequency  domain.  The  RML  provides  several  options  for  adaptive  delay 
estimation  in  the  time  domain.  In  one  configuration  the  RML  algorithm  Is  used 
as  an  adaptive  whitening  filter.  The  outputs  of  two  sensors  are  whitened  and 
then  cross-correlated.  In  another  configuration,  the  RML  algorithm  Is  used  to 
fit  an  MA  model  relating  the  output  of  one  sensor  to  the  output  of  another. 
This  Idea  Is  similar  to  the  one  proposed  in  [16].  Both  of  these  techniques 
seem  to  work  quite  well  as  reported  in  [1],  However,  a  more  interesting 
possibility  Is  to  use  a  multichannel  version  of  the  RML,  as  will  be  discussed 
in  section  3. 
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2.3  PERFORMANCE  EVALUATION 

The  single  channel  RML  algorithm  was  tested  in  the  context  of  several 
specific  applications,  as  described  above.  Based  on  our  experience  so  far  the 
RM.  appears  to  be  a  viable  technique  for  adaptive  HR  filtering  capable  of 
operating  at  low  SNR's.  The  RML  algorithm  provides  asynptotical ly  efficient 
parameter  estimates,  i.e. ,  the  variance  of  the  parameter  estimate  $  approaches 
the  Cramer-Rao  lower  bound  [17].  The  asymptotic  performance  bound  of  the 
algorithm  can  be  evaluated  using  the  results  presented  in  [18].  The  algorithm 
is  guaranteed  to  converge  when  the  signal  yt  is  an  ARMA  process,  provided  that 
the  orders  NA  and  NC  are  at  least  as  high  as  the  true  model  order  [19], 

[20].  The  SNR  gain  of  the  IIR-ALE  was  evaluated  analytically  in  [2].  The 
performance  of  the  RM.  in  other  applications  was  evaluated  mainly  by 
simulation  studies  [3],  [4],  [21]  in  which  it  compared  favorably  with 
currently  used  adaptive  1IR  filters  and  AR  estimation  techniques  such  as  the 
Maximum  Entropy  Method. 

The  main  difficulty  encountered  with  the  RML  algorithm  is  its  fairly  long 
transient  phase  which  typically  lasts  for  several  hundred  data  points.  In 
signal  processing  applications  involving  continuous  filtering  of  data,  the 
transient  phase  is  of  little  importance.  Once  the  algorithm  "locks-on"  the 
right  parameter  estimates  it  will  continue  tracking.  In  applications 
involving  short  data  records  the  convergence  rate  of  the  RML  may  not  be 
sufficiently  fast.  In  general  we  found  that  the  convergence  of  the  RML  was 
accelerated  as  the  parameter  k  in  the  modified  pre-filter  was  made  smaller. 
This  is  in  agreement  with  observations  by  other  researchers  that  the  Extended 
Least-Squares  technique  (k*0)  Is  usually  faster  than  the  RML  algorithm 
(k=l).  It  is  possible  that  further  modifications  can  be  made  to  increase  the 
convergence  rate  of  the  algorithm.  Unfortunately,  no  analytical  results  are 
available  on  the  convergence  rate  of  recursive  stochastic  algorithms  of  this 
type.  It  will  be  necessary  to  develop  appropriate  analysis  tools  before  this 
issue  can  be  properly  addressed.  However,  the  performance  of  the  RI4.  in  its 
pesent  form  is  probably  good  enough  for  most  adaptive  processing  applications. 

Finally  we  note  that  the  RML  algorithm  was  tested  mainly  on  stationary 
data.  Its  capability  for  tracking  signals  with  time-varying  statistics  was 
not  evaluated  adequately  due  to  time  and  budget  limitations. 
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3.  THE  MULTITARGET  ALGORITHM 


The  MTS  concept  Is  based  on  fitting  a  multi-input  multi-output  parametric 
model  to  the  observed  data.  In  [1]  and  in  appendix  B  we  show  that  the 
following  transfer  function  matrix  arises  from  the  physical  model: 

H(z)  =  M(z)D(z)  _1  (12) 

where 

D ( z )  ■  diag  { 1/d^ (z) }  =  spectral  model 

-T.  . 

N(z)  =  [z  1J]  =  propagation  model 

(t. .  =  propagation  delay  from  target  i  to  'sensor  j) 

^  J 

The  outputs  of  the  sensors  Y(z)  can  be  written  as 

Y ( z)  -  H( z)U(z)  +  V ( z)  (13) 

where  U(z),  V(z)  are  independent  white-noise  processes.  Using  this  problem 
formulation  we  need  a  way  of  estimating  the  parameters  of  the  parametric  model 
H(z)  from  observations  of  the  sensor  data  sequence.  Unfortunately,  the 
parameters  of  the  Right  Matrix  Fraction  Description  (RMFD)  H(z)  =  N(z)D(z)_1 
are  hard  to  estimate.  It  is  much  easier  to  estimate  the  parameters  of  a  left 
Matrix  Fraction  Description  (LMFD)H(z)  *  A"l(z)  B(z).The  RMFD  arises  naturally 
from  the  physical  model  and  the  target  parameters  appear  in  it  directly. 
Furthermore,  the  RMFD  has  a  very  special  structure  (D(z)  a  diagonal  matrix  and 
N(z)  a  delay  matrix)  while  no  such  structure  is  apparent  in  the  LMFD. 

However,  because  of  the  much  greater  complexity  of  estimating  RMFD  parameters, 
we  have  concentrated  so  far  on  the  LMFD. 

In  section  3.1  we  present  an  extension  of  the  RML  algorithm  for 
estimating  the  parameters  of  a  LMFD  of  the  transfer  function  H(z).  In  section 
3.2  we  present  some  simulation  results  using  this  algorithm.  The  problem  of 
uniqueness  and  other  theoretical  Issues  are  discussed  in  appendix  B. 
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3.1  THE  MULTICHANNEL  RM. 


The  RM.  algorithm  described  in  section  2  can  be  extended  to  estimate  the 
parameters  of  multichannel  ARMA  models  of  the  form 


pin  nu 

't  '  -£ai  Vi  ♦£bi  Vi  *  vt 


where  y,v,  are  p  x  1  vectors  and  A^ ,  are  p  x  p  matrices,  p  being  the  number 
of  sensors.  The  components  of  v^  are  uncorrelated  white  noise  processes.  Eq. 
(14)  can  be  written  as 

yt  =  9T  \  +  vfc  (15) 


where 


0T  =  [Ap...  Ana,  Bp  ... ,  Bnb]  ,  P  X  (NA  +  NB)p  matrix 

\  *  ”y t-NA ’  vt-l*"*  ’  ^NA  +  X  1  vector 


or  it  can  be  decoupled  into  p  equations  of  the  type 


y{ s  4  +  4* 


J  -  1 . .  (16) 


where  9^  is  the  j-th  column  of  9.  For  each  of  these  equations  we  can  now 
write  down  the  single  channel  algorithm  described  earlier.  All  of  these 
algorithms  will  have  a  common  gain  vector  Let 

*t  *  ^"yt-l’**’‘  -yt-NA*  ^-l*****  et-NB^  *17a* 

\  *  C-yt-l”*-*  ~yt-NA ’  ®t-l,‘”’  ®t-NB ^  *17b* 

where  et  yt,  et  are  to  be  defined.  The  update  equations  are  then  given  by 

pt  ■  [pt-i  -  Vi  *t  4  Vi/Ut  *  ♦»  \  ,18) 

«i  ■  H-i  *  pt  *t  <*t  •  *t  H-i>  1191 
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et  *  yt  ~  *t  H  •  et =  et-'T 

(20) 

yt  =  o^(z)  yt  ,  et  =  D"1  (z)  et 

(21) 

Dt(z)  =  1  +  kB1(t)z“1+  ...  +  KNBBNB(t)z'NB 

(22) 

The  RM.  algorighm  described  above  estimates  the  coefficients  of  the  prediction 
filter 

NA  MB 

•t  ’  H  n-f£8’t  “-11  H-i  1231 

i=l  i =1 

where 


yf  -  <>*9: 


tat- 1 


The  prediction  error  sequence  is  white,  but  there  is  no  guarantee  that  the 
components  are  uncorrelated.  Since  the  original  model  (14)  involved  a 

vector  with  uncorrelated  components,  we  need  to  decorrelate  the  components  of 
c.  as  well.  This  can  be  done  by  pre-multi plying  with  R”e/^,  where  R^2  is 
the  upper  triangular  square  root  of  the  prediction  error  covariance  matrix. 

In  other  words,  the  estimated  model  will  be 

H(z)  *  A^U)  Bt(z)  R^e/2  (24) 

The  square-root  covariance  matrix  R^  can  be  computed  recursively  using 
the  procedure  described  next.  First  note  the  following  fact: 


R^  =  ( 1  )  R^_j  +  w^  a  eQe^  (25a) 

wt  1  Vl/(wt-l  +  V  ’  wo  "  1  (25b) 

The  parameter  wt  is  needed  to  make  R*  an  unbiased  estimator  of  Eie^el).  If 

^  t  t  t 

we  were  to  compute 
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(26) 


R 


e 


h*L 


+  e 


T 

tet 


then  R*  will  be  a  scaled  version  of  the  true  covariance,  where  the  scaling 
factor  is  n 


Vt-1  +  1 


(27) 


For  example, if  \  ~  1  then  -n  =  t  and  the  unbiased  covariance  estimate  is 

i  W  v 

i  r^.  The  unbiased  covariance  estimate  R^/nt  obeys  the  following  recursion: 


^t  .  \  Vl  Rt-1  +  1_  T 
\  ”  \  \_i  cth 


3  (l-l/nt) 


(28) 


Let  wt  3  l/nt,  and  rename  R^/nt  as  R^  to  get  Eq.  (25a).  Note  also  that  Eq. 
(27)  can  be  rewritten  as  (25b). 


Instead  of  using  Eq.  (25)  and  factoring  it  at  each  time-step  to  get  the 
square-root  R^2,  it  is  possible  to  update  directly  R^2.  Note  that 


d£/2 

1  *t  Kt-1 

a  e/2 

) -  T 

/wt  et 

0 

J 

pre-array  post-array 


(29) 


where  ~  denotes  equi valance  up  to  an  orthogonal  transformation.  Note  that  the 
right-hand- side  matrix  Is  upper  triangular  and  can  be  obtained  from  the  left- 
hand-side  matrix  by  applying  a  Householder  trlangularlzation  routine  (see 
appendix  0  for  details).  The  complete  alorlthm  for  computing  R^2  is  as 
follows: 
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update  wt  (25b) 

~  i>\ 

conpute  /l-w  R^,  ^wt  et  and  set  them  in  a  pre-array 
apply  a  Householder  triangularization  routine  (L=p+1,  M=N=p) 
the  upper  p  x  p  matrix  in  the  post-array  is  R^2 

To  ensure  the  convergence  of  the  RML  algorithm  it  is  necessary  to  monitor 
the  stability  of  the  polynomial  det  B(z)  as  explained  in  section  2.1.  Target 
parameters  can  finally  be  obtained  from  the  estimated  transfer  function  matrix 
H(z)  (24),  as  explained  in  Appendix  3.  In  our  simulations  we  estimated  the 
spectrum  of  Target  i  by  computing  the  power  spectrum  of  the  impulse  response 
of  H(z)  from  its  i-th  input.  The  TDOA'  s  are  computed  by  cross-correlating 
impulse  responses  at  different  outputs. 

3.2  SIMULATION  RESULTS 

The  multi  target  tracking  concept  was  tested  by  computer  simulation.  Some 
preliminary  results  are  presented  here.  All  of  the  tests  so  far  were  for 
fixed  (non-moving)  targets  with  stationary  spectra.  The  RM-  algorithm 
described  in  section  3.1  was  used  to  estimate  the  ARMA  coefficients.  At  this 
stage,  no  attempt  was  made  to  enforce  the  special  structure  of  the  transfer 
function  H(z).  Therefore,  uniqueness  is  guaranteed  only  for  test  cases 
corresponding  to  stable  and  minimum  phase  systems.  As  discussed  in  appendix 
B,  the  more  general  case  requires  the  enforcement  of  the  special  structure  of 
N(z)  and  D(z). 

All  of  the  following  test  cases  were  run  on  N=1024  data  points  and  an 
input  signal  to  noise  ratio  (at  each  receiver)  of  20dB. 
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Case  1 


Target  Mo.  1 

Spectral  parameters:  dj(z)  =  1  -  1.64Z'1  +  Q.95z"2 
Location  parameters:  TDOA  =  2 


Target  No.  2 

Spectral  parameters:  d2(z)  =  1  -  1.34z-1,  +  0.95z“2 
Location  parameters:  TDOA  *  -2 


'1  .9z"2' 

V - 

O 

M 

»  ♦ 

_.9z~2  1 

.  0  d2(zj 

Fig.  4  depicts  the  cross-correlation  and  spectra  of  the  impulse  responses 
of  the  estimated  transfer  function  HU)  •  Note  tnat  the  correct  delay 
estimates  were  obtained  (the  peaks  of  the  correlation  function).  Note  also 
the  high  degree  of  spectral  separation  that  was  obtained!  The  estimated 
spectrum  of  target  No.  1  still  has  a  component  of  the  second  target,  but  that 
component  is  attenuated  by  more  than  30d8.  The  labeling  of  the  targets  (No.  1 
and  No.  2)  is,  of  course,  arbitrary  but  the  algorithm  correctly  associated  a 
spectrum  with  a  location  (TDOA). 


Target  No.  1 

Spectral  parameters:  dj^z)  =  1  -  1.64z‘*  +  0.95z"2 
Location  parameters:  TDOA  =  -2 


Target  No.  2 

Spectral  parameters:  d2(z)  *  1  -  0.95z-1  +  Q.95z-2 
Location  parameters:  TDOA  =  -2 


'l  .9z’2" 

d1(z)  0- 

,.9z'2  1 

.0  d2  ( z ) 

The  results  are  depicted  in  Fig.  5. 

Case  3 

Target  No.  1 

Spectral  parameters:  djfz)  =  1  -  1.39z_1  +  0.8z"2 
Location  parameters:  TDOA  *  2 


Target  No.  2 

Spectral  parameters:  d2(z)  *  1  +  G.8z-2 
Location  parameters:  TDOA  =  -2 


H(z)  * 


ri 


dl(z) 


0 

d2(z). 


-1 


These  targets  had  broader  spectra  then  the  targets  in  Case  1  and  Case  2.  The 
results  are  depicted  in  Fig.  6. 


2.00 


4.00 


Case  4 


Target  No.l 

Spectral  parameters:  ( z)  =  1  -  1.64Z"1  +  Q.95z'2 

Location  parameters:  TDOA  =  0 


Target  No. 2 

Spectral  parameters:  d2(z)  *  1  -  i.i4z“1  +  0.95z“2 


Location  parameters: 


H(z) 


1  .9z’2 
.1  .9z_1 


TDOA  *  -1 
d^z)  0 
0  d2(z). 


The  results  are  depicted  in  Fig.  7. 


In  all  the  test 'cases  so  far  we  have  chosen  N(z)  so  that  H(z)  is  a 
minimum  phase  plant.  When  H(z)  is  non-minimum  phase,  the  estimated  transfer 
function  H(z)  is  a  transformed  version  of  the  true  transfer  function,  i.e., 
H(z)  *  H(z)T(z),  T(z)  3  a  para-unitary  matrix.  Thus,  the  spectral  and 
location  parameters  can  not  be  read  directly  from  the  impulse  response.  The 
following  test  illustrates  this  fact. 

Case  5 

Target  No.l 

Spectral  parameters:  d^(z)  *  I  -  I.64z_1  +  .95z-2 
Location  parameters:  TDOA  *  1 


Target  No.  2 

Spectral  parameters:  d2(*)  =  1  -  1.34Z'1  +  .95z"2 
Location  paramets:  TOOA  »2 


The  results  are  depicted  in  Fig.  8.  Note  the  loss  of  spectral  seperation  and 
the  somewhat  ambiguous  delay  estimates. 
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-  Target  So.  1 

-  Target  So.  2 

Figure  7b.  The  es tinned  sDectra  of  two  targets 
N» 1024  data  ooints,  k*.9  . 


2. 


0.00  2.00  4.00 

-  Target  No.  1 

.  Target  No.  2 

Figure  3a.  Cross-correlation  of  the  imulse 

responses  corresponding  to  each  target 
N=1024  data  points,  k=1.0  . 


150.  290.  259. 


— -  Target  No.  t 

“  Target  No.  2 

Figure  3b.  The  estimated  spectra  of  two  targets 
N* 1024  data  points,  ks1.0  . 


The  results  presented  above  are  by  no  means  conclusive.  Only  a  limited 
number  of  test  cases  were  run  since  the  multichannel  program  is  difficult  to 
use  in  its  present  form.  Also,  further  work  is  needed  to  solve  adequately  the 

*  i  ~  -p/2 

uniqueness  problem.  The  LMFD  representation  At  (z)  Bjz)  Rt  '  has  to  be 
converted  into  an  RMFD  with  the  appropriate  structure  in  order  to  get  the 
correct  TOOA's  and  spectral  estimates  in  all  situations  (cf -  test  case  no. 

5).  Better  yet,  an  algorithm  for  estimating  directly  the  RMFD  needs  to  be 
derived. 
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4.  WORK  IN  PROGRESS 


In  phase  II  we  have  completed  the  development  of  the  single  channel 
algorithm.  While  some  Issues  such  as  convergence  rate  and  tracking  capability 
are  still  not  completely  resolved,  the  algorithm  is  basically  realty  for 
testing  on  real  data  and  for  a  conclusive  performance  evaluation.  This, 
however,  is  beyond  the  scope  of  the  MTS  project.  The  single  channel  algorithm 
will  not  be  pursued  further  at  this  point. 

The  main  objective  of  our  current  work  and  its  continuation  in  phase  III 
is  the  development  and  testing  of  the  multi  target  algorithm.  'While 
significant  progress  was  made  in  phase  II,  the  multi  target  algorithm  is  not 
yet  fully  developed,  as  was  discussed  in  Section  3.  We  are  currently  pursuing 
several  lines  of  investigation: 

(i)  Analysis  of  the  multi  target  case. 

The  nonuniqueness  problem  and  its  physical  interpretation  are  being 
studied.  We  are  looking  for  ways  of  estimating  correctly  target  spectra  and 
TDOA’s  from  the  LMFD  parameters.  The  problem  of  modeling  Doppler  effects  and 
incorporating  them  in  the  MTS  processor  is  being  considered.  Another  issue  is 
the  determination  of  the  number  of  targets  and  handling  non-square  transfer 
function  matrices. 

{ i  i )  Refinement  of  the  multichannel  Rfi.  algorithm. 

The  algorithm  described  in  Section  3.1  has  been  tested  under  high  SNR 
conditions.  Its  performance  is  severely  degraded  at  low  SNR’s.  Various 
methods  for  making  the  algorithm  more  robust  will  be  explored.  This  will 
require  the  enforcement  of  some  structural  constraints  to  reduce  the  number  of 
parameters  involved. 

(ill)  Alternative  methods. 

The  "ideal1’  approach  to  the  multi  target  problem  would  be  to  estimate 
directly  the  parameters  of  a  physical  model  (l.e.  the  RMFD,  for  the  situation 
assuned  in  this  report,  or  a  more  general  model  in  case  doppler  effects  are 
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included).  A  possible  approach  which  circumvents  some  of  the  difficulties 
involved  in  estimating  the  RMFD  parameters  is  to  work  with  the  spectral  matrix 
S(z)  *  H(z)  hT(z-I).  Another  possibility  is  to  develop  a  suboptimal  maximum 
likelihood  type  algorithm  for  estimating  directly  the  RMFD  parameters. 

Due  to  limited  resources  we  will  not  be  able  to  investigate  all  of  these 
Issues.  At  the  start  of  phase  III  the  most  promising  direction  will  be  chosen 
and  most  of  our  effort  will  be  devoted  to  it.  We  plan  to  develop  the 
multichannel  algorithm  sufficiently  so  that  by  the  end  of  phase  III  it  will  be 
possible  to  draw  conclusions  regarding  the  feasibility  of  our  parametric 
modeling  approach. 


5.  PROJECT  PUBLICATIONS 


This  section  lists  all  the  publications  resulting  from  the  Multitarget 
Tracking  Studies  project  to  date. 

Conference  Presentations  and  Publications 


1.  B.  Friedlander  and  J.J.  Anton,  “System  Identification  for  Multitarget 
Tracking",  Proceeding  of  the  13th  Asilomar  Conference  on  Circuits  Systems 
and  Computers,  Pacific  Grove,  California,  November  1979.  Also  presented 
to  the  National  Acadeny  of  Sciences  Panel  on  Applied  Mathematics  Research 
Alternatives  for  the  U.S.  Navy,  Washington  O.C.,  November  2,  1979. 

2.  B.  Friedlander,  "An  ARMA  Modeling  Approach  to  Multi  target  Tracking,"  Proc. 
of  the  19th  IEEE  Conference  on  Decision  and  Control,  December  1980. 

3.  B.  Friedlander,  "A  Pole-Zero  Lattice  Form  for  Adaptive  Line  Enhancement," 
Proc.  of  the  14th  Asilomar  Conference  on  Circuits,  Systems  and  Computers, 
Pacific  Grove,  California,  November  1980. 

4.  B.  Friedlander,  "System  Identification  Techniques  for  Adaptive  Noise 
Cancelling,"  Proc.  of  the  14th  Asilomar  Conference  on  Circuits,  Systems 
and  Computers,  Pacific  Grove,  California,  November  1980. 

5.  B.  Friedlander  and  M.  Morf,  "Least-Squares  Algorithms  for  Adaptive  Linear- 
Phase  Filtering,"  Proc.  of  the  1981  IEEE  Inti.  Conf.  on  Acoustics,  Speech 
and  Signal  Processing,  Atlanta,  Georgia,  March  30  -  April  1,  1981. 

6.  B.  Friedlander,  "A  Recursive  Maximum  Likelihood  Algorithm  for  ARMA  Line 
Enhancement,"  Proc.  of  the  1981  IEEE  Inti.  Conf.  on  Acoustics,  Speech  and 
Signal  Processing,  Atlanta,  Georgia,  March  30  -  April  1,  1981. 

7.  8.  Friedlander,  "A  Modified  Lattice  Algorithm  for  Deconvolving  Filtered 
Impulsive  Processes,"  Proc.  of  the  1981  IEEE  Inti.  Conf.  on  Accoustics, 
Speech  and  Signal  Processing,  Atlanta,  Georgia,  March  30  -  April  1,  1981. 


8.  B.  Friedlander,  "A  Recursive  Maximum  Likelihood  Algorithm  for  ARMA 
Spectral  Estimation,”  Proc.  Conference  on  Information  Sciences  and 
Systems,  Johns  Hopkins  University,  April,  1981. 

9.  B.  Friedlander,  "Appl ication  of  Recursive  Parameter  Estimation  Algoritms 
to  Adaptive  Signal  Processing,"  Proc.  Workshop  on  Adaptive  Processing, 

Yale  University,  Connecticut,  May,  1981. 

10.  B.  Friedlander,  "A  Modified  Pre-filter  for  Some  Recursive  Parameter 
Estimation  Algorithms”,  Proc.  20th  IEEE  Conf.  on  Decision  and  Control,  San 
Diego,  California,  December  16-18,  1981. 

Submitted  for  Journal  Publication 


1.  B.  Friedlander,  "Lattice  Implementations  of  the  Adaptive  Line  Enhancer," 
IEEE  Trans,  on  ASSP. 

2.  B.  Friedlander  and  M.  Morf,  “Least-Squares  Algorithms  for  Adaptive  Linear- 
Phase  Filtering",  IEEE  Trans,  on  ASSP, to  appear. 

3.  B.  Friedlander,  "A  Recursive  Maximum  Likelihood  Algorithm  for  ARMA  Line 
Enhancement,"  IEEE  Trans,  on  ASSP,  to  appear. 

4.  B.  Friedlander,  "System  Identification  Techniques  for  Adaptive  Noise 
Cancelling",  IEEE  Trans,  on  ASSP. 

5.  B.  Friedlander,  "A  Modified  Lattice  Algorithm  for  Deconvolving  Filtered 
Impulsive  Processes,"  IEEE  Trans,  on  ASSP. 

6.  B.  Friedlander,  "A  Modified  Pre-Filter  for  Some  Recursive  Parameter 
Estimation  Algorithms,"  IEEE  Trans,  on  AC. ,  to  appear  February  1982. 

7.  B.  Friedlander,  "System  Identification  Techniques  for  Adaptive  Signal 
Processing,"  IEEE  Trans,  on  ASSP. ,  to  appear.  Will  also  be  published  in 
the  Journal  of  Circuits,  Systems  and  Signal  Processing. 
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8.  B.  Friedlander,  "A  Recursive  Maximum  Likelihood  Algorithm  for  ARMA 
Spectral  Estimation,"  IEEE  Trans.  Information  Theory,  to  appear,  July 
1982. 

Reports 

B.  Friedlander,  "Multitarget  Tracking  Studies:  Phase  I  Final  Report," 
Report  No.  5334-01,  July  1980. 
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SYSTEM  IDENTIFICATION  TECHNIQUES  FOR  ADAPTIVE  SIGNAL  PROCESSING* 


B.  Friedlander 

Systems  Concrol,  Inc 
1301  Page  Mill  Road 
Palo  Alto,  CA  94304 


ABSTRACT 


Many  problems  in  adaptive  filtering  can  be  approached  from  the 
point  of  view  of  system  identification.  The  recursive  maximum 
likelihood  algorithm  is  proposed  for  estimating  the  parameters  of  the 
signal  model.  The  parameter  estimates  are  then  used  to  form  an  adaptive 
infinite  impulse  response  filter.  Several  examples  are  discussed 
including:  adaptive  line  enhancement,  adaptive  deconvolution,  adaptive 
noise  cancelling  and  adaptive  time  delay  esitmation. 


*This  work  was  supported  in  part  by  the  Office  of  Naval  Research, 
Contract  No.  N00014-79-C-0743. 


1.  INTRODUCTION 

Considerable  progress  -is  sad*  in  the  last  dacada 
in  the  development  and  analysis  of  racurtiva  parameter 
estimation  algorithms.  Tha  major  pare  of  this  work  was 
in  tha  araa  of  system  identification,  in  tha  context  of 
controlling  plant!  with  unknown  or  slowly  tiaa  varying 
paraaatars.  A  larga  ntnbar  of  algorithms  wara  davalopad 
for  fitting  linaar  nodal*  to  tha  observed  data.  Tha 
following  aRMAX  nodal  is  an  example  of  tha  class  of 
nodal*  chat  ara  typically  consldarad. 

Lac  u.  and  y.  danota  tha  input  and  output 
procassas  of  tha  nodil,  and  v.  an  unaaaiurabla  whlta 
nolsa  procass  (i.a.  a  "disturbance”,  in  tha  control 
terminology; .  Thasa  procassas  ara  ralatad  by  tha 
following  aquation: 

NA  MB  MC 

•  -  a  v  +■  b  «  ♦  ;  v  •  v  (2.) 

*  ’  ^  - 

which  can  be  -written  in  polynomial  fora  as 

A(:'‘iy.  •  3(c':iu.  *  C(*-l)v.  <Z*> 


3U*“)  -  b:a*‘  * 


s*1  *  unit  delay  operator,  i.a.,  t”ix,  •  x„^  .lei 

Tha  case  where  3(t  •  0  is  of  special  interest  in 

signal  processing  applications,  sine*  tha  output 
y.  *  (C(a“l)/A(**s'i  jv^  Is  than  an  autoragrassiva 

novipg-everage  (ARMAi  procass.  Whan  Stc*1)  •  3  and 
C(r”i')  *  0  -we  have  an  autoragrassiva  (AR)  procass 
y.  •  (1/A(*”’*))v  .  Such  procassas  ara  vary  coanon  in 

tlaa  series  analysis  and  statistical  signal  processing. 

Many  problems  in  signal  processing  involve  signals 
that  can  ba  represented  by  linaar  model*  of  this  c yp« 

(AR,  ARMA,  ARMAX),as  wa  will  show  later.  Knowledge  of 
tha  signal  nodal  paraaatars  aakes  it  relatively  scraignt- 
forward  to  design  filters  that  perform  various  process¬ 
ing  functions  such  as:  linaar  prediction/smoothlng. 
deconvolution,  noise/ Interference  suppression,  or 
spectral  analysis.  When  the  signal  nodal  is  not  xnown. 
the  parameters  of  tha  ralatad  filter  need  to  be  adap¬ 
tively  adjusted.  A  technique  that  is  cotoonly  used  in 
adaptive  control  problems,  is  to  first  estimate  tha 
nodal  parameters  end  than  design  tha  controller  as  if 
these  estimates  ware  the  :rue  parameter  values. 


The  same  ;au  can  c*  used  for  adaptive  signal  process¬ 
ing,  aa  airaadv  indicated  is  !1!-".-'. .  Th*  combination 
a  recursive  paraaacar  estimation  aigoricsm  with  a 
filtering  tecnniqu*  besec  an  known  modal  parameters  is 
the  aaih  cnasa  or  this  paper. 

The  successful  application  of  syscam  identification 
techniques  in  adaptive  control  motivated  us  to  apply 
the  same  techniques  to  signal  processing.  Of  particu¬ 
lar  interest  is  the  fact  that  some  commonly  used  param¬ 
eter  estimation  algoric.-a*  such  as  rtcursive  maximum 
livelihood  (S6C.1 .  extended  least-squares  (ELS)  or 
recursive  generalized  least  squares  (SOLS)  i  5] -[b] ,  are 
capable  of  estimating  ARMA  {or  ARMAR)  parameters,  and 
toe  just  AS  parameters,  as  ua  shall  sea,  this  leads 
naturally  to  adaptive  infinite  impulse  response  (HR) 
filters.  Adaptive  1IR  filtering  is  considered  by  th* 
signal  processing  community  to  be  a  difficult  problem. 
Consequently,  the  overwhelming  majority  of  the  work  in 
the  area  of  adaptive  signal  processing  seems  to  con¬ 
centrate  or.  finite  impulse  response  .'HR)  filtering, 
the  application  of  system  testification  algorithms 
opens  th*  way  to  th*  development  of  a  whole  new  class 
of  adapciv*  HR  filters,  banned  oy  the  extensive  con¬ 
vergence  analysis  that  was  performed  in  recent  years 


In  spit*  of  th*  natural  interconnections  between 
system  identification  and  adaptive  signal  processing 
very  little  work  semns  to  have  been  don*  to  transfer 
algorithms  from  on*  area  to  th*  ocher.  Th*  objective 
of  this  paper  is  to  report  on  some  recast  work  in  wcucr. 
a  version  of  the  RML  algorithm  was  successfully 
applied  to  a  numoer  of  problems  including:  adaptive 
line  enhancement, aaeptiv*  deconvolution,  adaptive  noise 
cancelling  and  time-delay  estimation.  A  brief  descrip¬ 
tion  of  these  applications  is  presented  in  Section  i. 

tn  Section  1  w*  present  th*  RJC.  algorithm  and 
discuss  its  properties.  Some  of  the  special  character¬ 
istics  of  the  signal  processing  problma  (when  compared 
co  th*  control  problem;  and  their  effect  on  the  behavior 
of  the  algorithm  ere  also  discussed.  Finally,  in 
Section  -  w*  outline  soma  alternative  algorithms  for 
adaptive  processing  and  areas  for  further  investigation. 
«•  hop*  tnat  auc  work  will  stimulate  further  research 
into  th*  numerous  potential  applications  of  parametar 
estimation  algorithms  to  adaptive  signal  processing. 

:.  TEE  RECURSIVE  MASSKISi  LIKELIHOOD  ALGORITHM 


Th*  RFC.  algorithm  provides  a  recursive  estimate 
of  th*  parameters  of  th*  ASMAR  model  in  equation  (1). 

For  a  detailed  derivation  of  this  elgoritia  tee  [13! 
[Hi.  Her*  we  present  only  a  jumaary  of  the  recursions. 


Let  i  denote  the  vector  of  model  parameters 

r  •  ci . esc’  > 

sod  ..  denote  th*  date  vector  and  the  filter  date 
vector,  respectively. 

•*.  -  ‘~?t-l . ’yt-SA’ut.l,'-*,u:-SJ’*t-l,'“,*t-NC11 

Ub) 

t  *  . ■>’t-KA,ut-l . “c-KI’*c-l,‘”,*c-HC'  * 


(3c> 

Denote  by  r,  the  parametar  estimates  sc  time  t,  end 
by  C(t  A)  th*  filter  whose  coefflcsnts  sr*  th*  estlmatei 
t, .  Then, 

*  V;<-1  -  *  prediction  error  (»*) 


V  m  '  9  ?  (  -  *  3 

•  c  '  c-1*  v  cr:-r] 

»  srror  covariance  aacri:*: 

•  ~  ^  p  ; 


y.  *  il/C(k*“l)>y.  } 
u.  -  (i/C(kz"i'-)u„  v  fi 


leered  quantities 


1/Ctka 


with  initial  conditions 

?  •  il,  i  »  initial  estimate  of  th*  covariance 


-  •  r  ,  initial  estimate  of  i  :  typical! v 

o  o 

-  •  0) . 


Th*  parameter  .  is  ad  exponential  weighting  on  th* 
date.  Typically  .  is  a  constant  close  co  unity,  or 


,.j  -  ,  •  0.M)  -5-i 

The  paremecer  k  is  used  to  "pull  in"  th*  roots  of  the 

polynomial  C\its”k)  into  the  unit  tirtle,  wnen  Hi''; 
hes  roots  near  the  unit  circle,  as  is  often  th*  case  is 
signal  processing  applications.  This  parameter  effects 
the  convergence  race  of  the  algorithm  as  discussed  in 
!  15!  •  la  fact,  for  k»0  this  algorithm  becomes  th*  so- 
celled  Extended  Least-Square  algorithm  detcribed  us  [ip[, 
)  1 ' ' .  In  most  cases  th*  choice  of  k  is  not  critical 
ana  in  ch*  following  we  assume  that  k  is  dose  or 
equal  to  unity.  To  ensure  convergence,  th*  stability 

of  CU*'-)  needs  tb  be  monitored.  If  unstable,  th* 
parameter  estimates  need  to  be  projected  into  a  region 
of  stability  [T.,*®;. 


Th*  asymptotic  properties  of  the  RML  algorithm 
have  been  investigated  in  considerable  detail.  It  was 
shown  chat  asymptotically  the  recursive  maximum  likeli¬ 
hood  technique  has  the  tame  properties  as  the  corres¬ 
ponding  off-line  version.  Thus  nothing  is  sacrificed 
by  going  to  t  recursive  implemencet ion .  provided  that 
enough  data  it  available.  Th*  maximum  likelihood  esti¬ 
mator  has  all  th*  desirable  properties  on*  may  expect 
from  t  parameter  estimator 

a  Aeympcotlc  consistency  il8!,  i.e.,  -  i 

as  th*  number  of  dace  points  V  goes  to 
inf inltv 

*  Asymptotic  efficiency  [1?[,  i.e..  th* 
estimation  error  covariance  approaches 
the  Cramer-Rao  lower  bound 

*  Th*  estimation  error  distribution  is 
asymptotically  normal  [19). 

The  convergence  properties  of  the  RML  algorithm 
were  studied  by  Ljung  ’,7],i81  end  others  £51 .  [Li ! .  It 
was  shown  that  this  elgoritia  will  always  converge  to  s 
local  maximum  of  th*  likelihood  function.  In  son*  situ¬ 
ations  there  may  exist  "falsa"  maxima  which  can  cause 
difficulties.  However,  for  ARfiA  processes  ic  we*  shove 
that  ell  th*  local  maxima  coincide  with  the  global  maxi¬ 
mum  [DO!  (provided  that  the  orders  (Ka.NC)  of  the  esti¬ 
mated  A1HA  model  are  equal  to  br  larger  than  the  true 
model  orders). 


Relatively  little  is  known  about  ch*  convergent* 
race  of  th*  RXL  algorithm  for  different  type*  of  pro¬ 
cesses.  Hardly  any  analytical  units  are  available  end 
most  of  tb*  remiltt  are  based  on  extesalve  simulation 
studies  ill}.  However ,  most  of  that*  studies  were 


i.i 


celateo  t o  control  orobiem*  involving  stac-e  events 
vi:.-  poies  ana  ceroes  veil  inside  the  unit  circ.e. 

".odels  arcsine  in  scgnsl  processing  applications  f-pi- 
pall'-'  invoice  narrowoanc  signals,  wr.icr.  act  represented 
by  poles  and  : spots  pn  or  vary  near  tn*  unit  circle. 

3u r  oun  experience  indicates  that  pole  and  pare  loca¬ 
tions  have  a  significant  effect  on  convargance  pacts, 
i'han  :p. a  poias  and  raroas  ara  wall  saparatad,  fast 
ponvarganca  was  ooaarvad.  Clusters  of  solas  and  zeroes,  .&■ 
especial!'-’  - n a r  thay  ara  naar  tna  unit  circle,  ocean 
lead  to  much  slower  ponvarganca.  The  casus  of  con  var¬ 
iant*  pact  cs  crucial  in  adaptive  signal  processing 
since  signals  ara  often  nonatationary  and  tiae- 
varying  and  it  is  important  to  enow  how  wail  the 
adaptive  filter  will  trace  the  changing  parameters. 

In  adaptive  control  problems  the  parameters  ara  often 
varv  slovlv  tiaa  varying  .compared  to  the  tiaa  con¬ 
stants  of  the  plant;.  *a  ara  currently  investigating 
the  convergence  rata  of  the  RML  algorithm  tor  ARMA 
nodals  with  poles  and  leroea  naar  the  unit  circle. 

The  propartias  of  the  maximum  likelihood  estimator 
described  above  make  it  vary  attractive  for  adaptive 
IIS  filtering,  as  illuscraced  in  the  next  section. 


:  .  ADAPTIVE  SIGNAL  PROCESSING:  SOME  EXAMPLES 

3.1  the  Adaptive  line  Enhancer  ',AL£1 

the  ALE  is  an  adaptive  filter  for  narrowband 
signals  ir.  additive  noise  Ill’  •;  131.  The  .ALE  can  be 
interpreted  as  an  adaptive  predictor,  i.e.,  ita  output 
is  y.  .  ■  the  estimate  of  y  begad  on  dace  us  to 

tiaa  t-l.  A  narrowoand  'autoregrasslvai  signal  in 
wnita  noise  tan  be  represented  as  an  ARMA  proctss  IJij. 
Thus,  tna  optimal  predictor  la  given  by 


NA  SC 


which  is  cepicted  as  a  tapped  delay  line  filter  in 
Figure  1.  The  parameters  of  the  filter  will  be 
adjusted  using  the  RML  algorithm  ( with  ,VB»0).  Sot* 
that  th*  resulting  filter  has  an  inflnits  impulse  res¬ 
ponse  while  most  ALE-s  discussed  in  the  literature  are 
of  the  MR  type.  In  121]  we  discuss  in  dsteil  th# 
advantages  of  the  HR- ALE  and  its  superior  performance 
at  low  signal -to -no isa  ratios. 

To  illustrate  th*  behavior  of  the  IZS-ALE  va 
depict  in  Figure  1  the  input  and  output  of  the  filter 
for  a  single  sinusoid  In  noise,  at  SNR  •  0  dj.  The 
spectra  were  obtained  ay  a  312  point  FFT.  The  ALE 


is#; 


i 


Figure  lb  Spectrum  of  the  ALE  Output 

Z  Z  “A 

cutput  data  corresponds  to  the  last  311  samples  in  a 
total  of  1043  data  points.  Note  the  significant  noise 
reduction.  Further  noise  suppression  is  achieved  if 
the  filcer  is  allowed  to  continue  ira  convergence. 

5a*  12*1, 115)  for  more  raiulca. 

The  parameters  of  th*  AfttA  model  Kr’1)  ,'At.c"'". 
which  serve  as  th*  coefficient*  of  :n#  11S-A1E  can  also 
“**?  co  **£1“£*  £h#  spectrum  of  tn*  obssrved  signal. 
In  '  2b]  w*  present  some  comperlsons  of  th*  AJOiA  spectr* 
obtaineo  by  th*  KML  with  corresponding  estimates 
obtained  by  th*  mexvnmn  entropy  method. 


Adaptive  Deconvolution 


Th*  need  to  extract  *  signal,  given  a  filtered 
version  of  th*  signal  arises  in  many  situation*  incluo- 
ing:  (i)  speech  enelysis/synthesis  by  linear  predictive 
coding,  sad  pitch  estimation,  iii)  estimation  of  th* 


reflectivity  sequence  in  seismic  dace  processing,  till) 
channel  equal ice t ion  for  th#  removal  of  th#  interivsboi 
interference  caused  by  convolution  of  the  message 
sequence  with  th*  channel  impulse  response. 


Consider  th*  case  where  a  whit#  signal  process 
pastes  through  in  HR  filter  C(i<)/A(s**.'  Stable  and 
Biaimia  phase) .  Th*  SML  algorithm  can  be  used  to 
estimate  th#  filter  parameters,  end  deconvolution  will 
b*  achieved  by  peseingth*  data  through  th*  estimated 

inverse  filtsr  AU*l)/C(t*1) .  Not*  that  most  eurrsnt 
deconvolution  techniques  are  limited  to  th*  case  where 

th*  convolving  filter  has  only  poles  C./AU*1))  while 


cn*  KU.  can  handle  the  pole-sero  css*  (CU  1/AU  *)). 
In  Figure  3  w*  present  a  comparison  between  HR  and  n 
deconvolution.  The  data  y  was  generated  in  this  csss 
by  passing  *  train  of  impulses  through  th*  filter 
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Figure  1  Th*  HR  Adaptive  Lin*  Enhancer 


and  adding  measurement  noise.  Th*  signal  to  noise  retie 
v»#  20  d>.  Not*  that  the  HR  filter  restores  th* 


”5  mg  ch«  RML  algetic  bin  to  esciaace  :.n«  ARMAJ*;  aodai 
pa  rant  cars  anc  then  escinacing  toe  "clear."  signal  by 
aquations  '?>.  11'  gives  an  X1R-ANC .  Figure  -  depicts 

me  noise  cancellation  achieved  tor  a  narrowband  sicnai 
in  me  presence  of  narrowband  interference.  See  .21'. 
for  sore  results. 

c  a  .u  i44  u* 


•i.h  -  -  ■  —  --a.ee 

j«e.  3<s.  :*•.  «ee*  (Saaple  M«*ber; 

Figure  3b  The  Deconvolved  Signal  :  at  the 
Output  of  an  FIR  Filter  'VNA*2 

original  impulse  sequence  while  the  FIR  filter  gives 
oairs  of  impulses.  This  is  caused  by  ;he  fact  that  the 
moving  average  part  of  the  convolving  filter  is  not 
identified  by  the  algorithm  and  therefore  the  de- 

•1  ^ 

convolve  signal  it  Ait  )y.  •  Cf.i  >v.  -  t.I-z  ' ,v  . 

For  a  sort  dacallad  discussion  of  adaptive  daconvolu- 
:ion  t««  ;:7;,  ;:8;. 


Adapt  iva  Sotw  Cancelling 


Figure  «b  Spectrum  ot  ANC  Output  s 


Th*  ANC  and  its  applications  art  dlaouaaad  is  ;i9) 
.'30) .  Kara  va  pratanc  only  a  vary  briar  description. 

In  th*  ANC  problem  va  ara  providad  with  a  oolsy  aaas- 
uramanc  y.  ot  tha  signal  t, 

Tt  m  ».  *  *.  C3) 

and  alto  vith  a  reference  input  u,  vhich  contains  in¬ 
formation  about  tha  nolsa  process's..  Fro*  th*  tidt 

information  (u.)  an  estimate  z .  of  tha  noita  proeass 

is  abtainad  and  than  subtractad  iron  tha  primary  input 
?.  to  "cancal  out"  tha  aolsa,  i.a.. 


Not*  that  in  addition  to  obtaining  th*  pradlcco; 

5(r”1)/A(z'1')  of  th*  aois*  procast  z.,  tha  SKI  aljo- 
rltha  automatically  providas  ua  vith  in  AKXA  nodal 

for  th*  signal  proeass  s.  .  If  th* 
signal  s.  is  a  narrovband  process  eorruptad  by  vhita 
nolsa  (uncorralatad  with  th*  aolsa  procaas  i.  ju- 

turad  by  th*  reference  input) ,  additional  nolsa  sup¬ 
pression  can  b*  obtained  by  narrovband  filtering.  Thus, 
th*  AIWA  parameters  (C.A)  can  b*  used  to  form  an  IIR- 
ALE  at  was  discutsad  in  Section  3.1,  i.a. , 


under  th*  aasiapcions  that  u,  and  t  ara  ralatad 

<•  c 


by  a  linaar  modal  and  that 


is  an  aRMA  process  it 


can  b*  shown  [31J  that  7,  is  an  ASM*!  proeass  of  tha 

'am  * 


S(x  *) 

■  t.  •  ■'  _■  U.  -  ■ 
-  A(z  -)  ' 


*t  ‘  -  *1*1*  -  Vt-i  •  ;1:) 

i-i  1  i-1  1  E  1 

In  otharvords ,  th*  PKL  algorithm  lead*  naturally  to  an 
adaptive  filter  (aquation*  (9),  (11),  (12))  chat  can 
b*  interpreted  at  an  HR- ANC  followed  by  an  IIR-aLE. 

In  [31]  w*  praasmt  a  much  mot*  datallad  discussion  of 
this  filter. 

3 .1  Adaptive  Time  Pal*. 


vhar*  u.  is  th*  reference  Input  and  v  it  a  vhita 

-  C 

nolsa  procaas.  Tha  noise  attiaata  can  b*  obtained  by 
th*  following  III  filter. 


Tha  need  to  attiaata  cla*  delay  between  two 
signals  srisas  in  many  applications  such  as  target 
localisation  by  sonar  ty stoat,  and  position  estimation 
by  radio  navigation  systems.  Th*  problaa  it  uaually 
formulated  at  fallows:  lot  us  asanas  that  two  tensors 


.raeivers,  microphones,  geophones)  receive  time  shifted 
and  scales  versions  at  mm  signal  a.  *  *  n.  , 

y  •  dx.  -  a_  ,  where  a.,  a.  ar«  independent  mea¬ 
surement  noise  prec*s*«(.  “any  different  techniques 
far  estimating  the  delay  t  have  baas  propoaad  is  the 
literature.  These  technique*  typically  involva  filter¬ 
ing  the  signals  and  than  craas-carr slating  ]321-i34]. 

To  do  tills  in  as  aptiaal  way  (m  tha  leaec-equaraa  or 
saxinum  likelihood  sansa)  raqulras  knowledge  of  tha 
statistics  af  both  tha  signal  and  noiae.  '.’sing  tha 
system  identification  approach  we  wars  sola  to  derive 
an  adapclva  tacnaiqua  for  tins  dalay  astiaatlon  which 
raqulras  no  prior  information  [37], 133].  Hers  wa 
prasanc  only  tha  simplest  varsion  af  tha  algorltisas 
prasancad  in  137],  [35]. 


5quars-Sooc  Fora 

Tha  update  aquation  1 )  for  tha  arror  covariance 
natrix  ?.  suffers  from  numerical  problems  whan  the 

number  of  estimated  paraaatars  n  •  :<A  -  N'B  »  SC  is 
large  ie.g.,  a  ■  10).  A  such  batter  implementation  is 

obtained  by  using  tha  square  root  fora  in  which  ?*' * 

1/”  -/**  * 

is  propagated  rather  than  ?,  'where  ?.?*•?_. 

Tha  advantages  of  square-root  algorithms  wars  discussed 
in  daail  bv  31eraas  [37], 

"Test"  laplaaantation 

Tha  SML  algorltia  prasancad  in  aquation  (A)  ra- 


Soca  that  two  procasaas  that  are  delayed  versions 
af  tha  same  underlying  signal  are  related  by  a  aoving- 
average  filter  whose  coefficients  contain  tha  dalay 
inforaacion.  In  tha  example  mentioned  above 

y.  -  3(c*1)u.  *  v£  (.13) 

wnara  , 

S(c*i)  •  da"' 

Co,*Q  for  ii*T  and  b,-d  far  i«c) 

v  •  a  -  dn  •  a  white  noise  process 
•  .  t-T 

Using  tha  ML  algorithm  to  estimate  the  coefficients 

of  Sic"1"  will  provide  an  astiaata  of  tha  time  delay 
by  looking  at  the  index  af  tha  largest  coefficient  af 

.  If  tha  dalay  la  not  an  intagar  multiple  of 
the  sampling  incarval,  it  is  necasaary  to  perform  a 
simple  incarpolatlan  to  gat  at  tha  true  daisy  [37]. 

Tha  3 20.  algorithm  is.  of  course,  capable  of  handling 
the  caaa  where  v.  is  nor  white. 


A  much  sort  sophisticated  approach  is  based  on 
tha  idea  of  fitting  a  multichannel  modal  to  a  vector 
of  sansar  measurements  y,,  e.g.. 


SA  N'C 

7.  •  - 1  v— i 

i«l  *  e  1  i-1  *  '  1 

(11) 

wher* 

y.  is  a  p  *  1  vector  and  A, , 

Ci 

?  *  p 

matrices.  In  [27],  [36]  we  have 

show  that  the 

resulting  multi-input  multi-output  (xmo)  adaptive 
filter  can  be  interpreted  as  a  combination  of  an  adap¬ 
tive  baamforoar  and  a  time  dalay  estimator.  Tha 
interesting  point  brought  out  in  [36]  la  chat  this 
filtar  is  capable  of  handling  simultaneously  several 
(up  to  p;  targets.'  Sara  wa  only  note  that  tha  usa  of 
ittMO  signal  aodals  opens  tha  way  for  developing  naw 
classes  of  MXMO  adaptive  filters  which  have  numerous 
applications  in  array  processing,  beamforming,  and 
processing  of  multlaanaor  data.  Current  adaptive 
filtering  techniques  ara  almost  amtlrely  devoted  to 
single  input  single  output  filters.  We  believe  that 
parhape  tha  aoet  important  contribution  of  our  model- 
lag  approach  is  chat  it  provides  a  syscwacic  frame- 
wrk  for  handling  multichannel  problems. 

1.  CONCLUSIONS 

We  prasancad  an  approach  for  darvaloping  adaptive 
filters  for  various  signal  processing  problems .  While 
tha  ML  algorithm  described  in  this  papar  is  wall 
know,  its  application  to  tha  class  of  problama  de¬ 
scribed  is  Section  3  la  apparwtly  new.  The  ML  algo¬ 
rithm  wa  presented  bare  in  one  particular  form.  It 
la  important  to  noea  that  alternative  forma  can  ba 
uaad  to  derive  ocher  implements cions  at  chase  adaptive 
filters  with  similar  asymptotic  properties. 
examples : 


quires  in  the  order  of  »n*  -  So  multiplications  end 
additions  per  data  point.  For  higher  order  models 
(i.e.,  large  values  of  n)  the  amount  of  computation  may 
become  excessively  large.  Thus,  it  is  important  to 
search  for  mors  efficient  estimation  techniques,  ’-’sing 
tha  idea  of  "shift  low  rank"  developed  by  Morf  led  to 
Implementation!  of  the  ML  requiring  proportional  to  a 

rather  than  n"  multiplications  and  additions  [38]. 

This  technique  provides  an  efficient  wy  of  computing 
the  gain  vector  .  The  rest  of  the  algorithm 

remains  unchanged.  The  detailed  update  aquae  ions  can 
ba  found  in  [38]  and  will  not  be  repeated  here. 

Recursive  Lattice  Forms 

The  recently  developed  square-root  normalised  lat¬ 
tice  forts  [ 39]  combine  the  good  numerical  behavior  of 
square-root  implementation  with  computational  effi¬ 
ciency.  They  furthermore  provide  a  computational 
technique  chat  is  recursive  both  in  time  end  in  model 
order.  In  fact,  the  lattice  form  provides  simultaneous 
parameter  estimates  corresponding  co  filters  of  all 
orders  up  co  a  maximum  order!  This  is  very  sueful  for 
addressing  the  order  determination  problw,  which  is 
one  of  the  more  difficult  aspects  of  signal  modeling. 

Symmetric  Afa~S  Polynomial 

In  soma  situations  the  parameters  of  tha  ARMA  model 
sre  interrelated  in  soon  wy.  For  example,  the  case  af 
sinusoids  in  white  noise  can  be  show  co  have  a  sym¬ 
metric  A (a-1)  polynomial,  i.e.,  a  ■  e  In  ]*0]  we 

1  n  “V 

presented  e  wy  of  incorporating  this  constraint  in  the 
parameter  estimation  algorithm.  Several  high  resolution 
spectral  estimation  techniques  are  implicitly  "Sym¬ 
metrising"  the  predictor  coefficients.  In  general, 
whenever  the  problem  has  tome  special  structure  chac 
can  ba  used  to  reduce  tha  number  of  estimated  parameters, 
one  should  explore  the  possibility  of  using  that  struc¬ 
ture  in  Che  parameter  estimation  algorithm. 

Finally  we  should  note  chac  the  results  prasancad 
hara  ara  poly  preliminary.  Much  work  remains  to  be  done 
on  the  analysis  and  performance  evaluation  of  the  ML 
and  related  algorithms  in  adaptive  signal  processing 
applications.  Some  specific  issues  chac  need  to  be 
addressed  ere:  the  convergence  rate  of  the  ML  for 
4*4feraat  classes  of  signals,  tha  tradeoff  between 
Parameter  cracking  capability  (I.e.,  tha  value  of  ■) 
end  filtar  performance,  end  tha  development  of  robust 
techniques  for  ASIA  order  determination  in  real -time. 
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ABSTRACT 

A  new  approach  is  presented  for  locating  aultiple 
targets  from  signals  received  by  a  number  of  spatially 
distributed  sensors.  A  multi-input  multi-output  model 
is  assumed  to  generate  the  observed  data.  It  is  shown 
chat  the  model  parameters  are  related  to  che  locations 
of  the  targets  and  to  their  spectra.  The  structure  of 
che  model  is  explored  and  a  specific  algorithm  is  given 
for  estimating  its  parameters  from  sensor  measurements. 
The  sain  feature  of  the  proposed  approach  is  that  ic 
estimates  simultaneously  che  parameters  of  all  che 
targets  of  incerest. 

1.  INTRODUCTION 

Tracking  multiple  targets  represents  special  dif¬ 
ficulties  since  there  can  be  uncertainties  associated 
vlch  the  measurements  beyond  their  inaccuracy  which  is 
usually  modeled  by  some  additive  noise.  This  additional 
uncertainty  is  related  to  che  origin  of  the  measure¬ 
ments.  Since  several  targecs  are  present,  it  is  neces¬ 
sary  to  sort  ouc  which  measurement  corresponds  to  which 
target.  In  ocher  words,  in  addiclon  to  the  problem  of 
detection  and  bearing/ range  estimation,  there  is  a  pro¬ 
blem  of  properly  labeling  the  sec  of  measurements.  The 
latter  problem  is  usually  referred  to  as  target  associ¬ 
ation  or  crack  formation. 

Typically,  these  two  facets  of  multicarget  track¬ 
ing  are  treated  separately,  first,  a  sec  of  potential 
target  locations  is  obcained.  Then  some  method  is  used 
to  label  these  locations  by  the  targets  to  which  they 
correspond  in  a  manner  consistent  wlch  previous  measure¬ 
ments.  Techniques  for  labeling  or  multlcargec  tracking 
have  been  developed  using  various  approaches  including: 
Kalman  filtering  (for  active  sonar  ill;  for  radar  ill); 
dayesian  mechods  [31— [7];  Integer  programming  [31;  and 
Track  splitting  [9],  [10],  [1],  [I].  for  a  recent 
survey  see  [11J. 

In  all  of  these  techniques  the  basic  detection  and 
iocaclon  estimation  are  performed  separately  for  each 
target.  The  multicarget  aspect  of  che  problem  enters 
only  through  the  labeling  procedure.  In  ocher  words, 
the  tracking  problem  is  treated  as  a  collection  of 
single  cargcc  problems,  which  has  co  be  put  cogecher 
in  a  systematic  and  conslscanc  way. 

In  this  paper,  we  attempt  co  cackle  directly  the 
multichannel  nature  of  che  problem.  Instead  of  looking 
at  one  target  at  a  time,  we  want  to  estimate  simultan  - 
eouslv  multicarget  parameters  (such  as  location  and 
spectrum).  The  approach  Is  best  understood  in  the  con¬ 
text  of  che  passive  cracking,  where  an  array  of  sensors 
measures  signals  (electromagnetic,  acoustic  or  seismic) 
generated  by  targecs.  This  type  of  problem  arise*  in 
sonar  systems,  acouscic  surveillance  systems  (detection 
of  low  flying  aircraft)  and  seismic  arrays  (oil  explor¬ 
ation,  earthquake  localization.  Intruder  detection  or 
artillary  localization).  The  active  cracking  problem 
(radar,  active  sonar)  can  be  hendled  in  chis  framework 
but  will  not  be  discussed  here. 

The  proposed  spproech  is  bated  on  che  idea  of  fit¬ 
ting  a  multi-input  multi-output  model  to  the  vector 
time-series  observed  at  che  outputs  of  the  sensors. 

Under  certain  aeaumpclons  cha  parameters  of  this  mcdel 
can  be  shown  co  concain  Information  about  che  locations 
of  all  the  targets,  as  wall  as  other  useful  Information 
such  aj  target  spectra.  If  the  parameter  estimation  is 


performed  recursively  and  continually,  ic  is  hoped  that 
the  relationship  beeveen  che  model  parameters  and  the 
targets  co  which  <-hey  correspond  will  be  consistently 
maintained.  Once  che  targets  are  labeled  in  a  particu¬ 
lar  way,  this  labeling  will  stay  fixed  without  need  for 
rechecking  or  relabeling.  This  will  eliminate  the  need 
for  a  separate  step  of  target  assolcaclon  which  is 
inherent  in  current  multitarget  cracking  techniques. 

In  Section  2,  we  present  the  basic  problem  formulation 
and  explore  the  relationship  beeveen  model  parzmecers 
and  targec  location. 

The  formulation  of  target  tracking  as  a  multichannel 
signal  modeling  problem  raises  a  number  of  interesting 
questions  In  the  areas  of  3yscem  identification  and 
structure  of  multivariable  linear  systems.  Some  of  these 
issues  are  explored  in  Section  3,  and  some  in  Section  5. 
In  Section  4,  we  present  a  specific  algorithm  for  per¬ 
forming  parameter  estimation.  This  algorithm  was  used 
in  s  preliminary  simulation  study,  which  is  described 
in  che  Appendix.  This  algorithm  is  presented  here  only 
as  an  example,  and  should  not  be  interpreted  as  the  best 
choice  for  this  application. 

Finally,  we  note  that  the  emphasis  in  chis  paper  Is 
on  developing  a  new  framework  for  handling  the  multi- 
target  cracking  problem  and  exploring  some  of  the  theor¬ 
etical  issues  it  raises.  No  claims  are  made  regarding 
che  relative  merits  of  this  approach  to  current  track¬ 
ing  algorithms.  At  this  stage,  we  are  mainly  interested 
in  creating  directly  the  multichannel  aspects  of  che 
tracking  problem,  and  not  in  making  a  comparative  evalu¬ 
ation  with  current  techniques. 

2.  THE  MODELING  FRAMEWORK 

To  illustrate  the  basic  ideas  of  our  approach  we 
consider  first  the  single  target  case  depicted  in  Figure 
1.  Two  sensors  are  measuring  the  signal  x(0  propagat¬ 
ing  from  a  target  located  somewhere  In  che  plane.  V* 
assume  that  the  propagation  involves  only  some  time 


delays  and  attenuation.  Thus,  the  outputs  y^,  y, 
the  two  sensors  can  be  modeled  as 

of 

File) 

•  x(t  -  t2)  +  v^fc) 

(la) 

y,(t> 

-  nx(t  -  T,)  +■  v2(c)  , 

(lb) 

where  T^, 

t2  are  the  propagation  delays  from  che 

car- 

get  to  the  two  sensors,  n  represents  attenuedon  and 
v. ,  v,  are  independent  measurement  noise  processes. 

The  time  sampled  version  of  these  oueputs  will  be 
written  as: 

y^k)  -  x(k  -  Tx)  +  vL(k)  (2s) 

y2(k)  •  nx(k  -  T;)  +  v2(k)  ,  (2b) 

where 

t  •  kAT,  r.  •  TjAT,  T,  -  TjAT. 

Not*  that  the  delays  7^,  %  art  assumed  to  be  integer 
multiples  of  ths  sampling  period.  The  ess*  of  non- 
integer  delays  is  discussed  later. 
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Figure  1.  Two  Sensors  nod  One  Target 

In  many  applications,  the  singal  x(c)  can  be 
adequately  modeled  as  an  autoregressive  (AS)  process 
of  che  fora 


x(k-i)  +  u(k) 


where  u(k)  is  a  white  noise  process.  This  Is  partic¬ 
ularly  crue  for  narrowband  signals.  The  no re  general 
case  of  signals  with  rational  (ARMA)  spectra  is  dis¬ 
cussed  in  Section  5.  Taking  Che  z-cranafora  of  £q.  (0) 
we  gee 

y,(z)-z*rlX(z)+V1(z)-z'Tl'a/d(z))U(z)-t-V1(z)  (4a) 

v,(z)-nz  ^-X(z)+V,(z)”nz  *l(l/d(z))U(z)  +  V,(z) 

(4b) 

where  d(z)«l+d,  z  1  +■  ...  +  d?z  p.  Written  in  vector 
fora  the  transfer  function  from  u  to  Che  sensor  out¬ 
puts  y,,y,  is  given  by  M(z)/d(z),  where 

y(z)-N(z)X(z)+V(z)-N(z)(l,'d(z))b’(z)+V(z)  (6a) 


or  using  the  notation  y(iAT)  ■  y(i) . 

» 

y(i)  *  2  n(k)  x(i-k)  +  v(k) 

k*-» 

where  n(k)*sinc[(2-k)AT+df] .  The  coefficient  n(k) 
achieves  its  largest  value  for  k*l  or  111.  As  an 
approximation  to  (9;  we  may  consider  replaicng  Che 
infinite  sum  by  a  finite  sum 


y(i)  *  Jr  a(k)  x(i-k)  *■  v( 


Thue,  Che  polynomials  S^(z)  and  ’H^z)  in  (6)  have 
to  be  replaced  by 

a  . 

N,  (z)  -  £  n.(k)z  .  i  •  1,  2. 
k-0 

The  delays  T^.T^  can  still  be  estimated  by  looking  ac 

the  index  of  the  largest  coefficient,  or  by  using  inter¬ 
polation.  A  reasonable  estimate  of  Che  delay  can  be 
given  by  the  value  '  which  maximizes  the  function  n( ' ) 


k)  sine  (T-k) 


-..'oca  that  the  numerator  of  chis  transfer  function  con¬ 
tains  information  about  the  relative  delays  from  the 
target  to  the  sensors  (which  is  directly  related  to 
target  location) .  while  the  denominator  concains  spec¬ 
tral  information.  Thus,  if  we  could  estimate  the  para¬ 
meters  of  the  model  X(z)/d(z)  from  che  observed  data 
\y,(t),  y,(t)>,  we  could  easily  obtain  che  information 
needed  for  tracking.  As  long  as  we  do  not  have  direct 
measurements  of  x(t)  or  u(c)  the  only  information 
that  can  be  extracted  from  the  spectral  factorization 
problem  3y(i)*(X(z) /d(z) )  OUz'^/diz*1) )  ’  is,  of  course 
nonunique.  Arbitrary  delays  z”  can  be  inserted  into 
the  spectral  factor  wlchouc  changing  the  spectrum. 

Thus,  as  intuitively  expected,  the  absolute  delays 
Tj_,T,  cannot  be  determined,  however,  che  difference 
in  the  degreee  if  the  polynoaiials  M,(z),  X„(z)  will 
be  unique,  and  provide  information  aoout  the  time- 
dlffarenca  of  arrival  (TDOA) ,  Aj^T^-Tj.  TDOA's 

for  two  or  more  pairs  of  sensors  will  uniquely  locate 
the  target  [12]. 

The  case  where  propagation  dalaya  (or  chair  dif¬ 
ferences)  are  non-integer  multiples  of  the  stapling 
period  can  be  approximated  by  an  aucoregreaalve  moving- 
average  (ARMA)  model  similar  to  equacloo  (6) .  To  see 
chis  consider  the  signal  x(c+T),  and  a  sampling  inter¬ 
val  AT.  Then. 

x(c-n)  ■  T  x(kAT)  sine  (c-t'T  -  kAT)  (7a) 


slnc(t)  -  sin  UIc/AD/CIt/AT)  .  ( 

Let  T  ■  IAT  +  At  ,  0  <  Ar<AT,  and  let 

y(t)  •  x(t+t)  +  v(t)  .  ( 

Then  the  sampled  version  of  y(c)  will  be  given  by 


The  discussion  above  requires,  of  course,  thac  the 
sampling  race  1/AT  be  sufficiently  large  in  relation  to 
che  signal  bandwidth.  The  order  of  m  of  the  moving 
average  (MA)  part  of  the  model  N(z)  will  be  determined 
by  the  sampling  interval  AT  and  by  the  expected  spread 
of  TTJOA's.  The  latter  depends  on  the  maximum  change  in 
target  bearing  expected  over  che  processing  interval. 

It  is  possible  to  reduce  this  order  by  removing  "bulk” 
delays  prior  to  che  parameter  estimation  algorithm. 

The  order  p  of  the  AR  part  of  the  model  d(z)  will 
depend  on  the  number  of  spectral  lines  in  che  target 
spectrum  (e.g.,  twice  the  number  of  important  lines). 

In  order  to  reduce  the  required  sampling  rate  1/AT 
and  the  modal  order,  it  is  assumed  che  signals  have 
been  shifted  to  a  baseband  (e.g.,  removal  of  carrlar. 

If  pr8sanc) .  For  a  more  detailed  description  of  9oma 
of  the  practical  Issues  involved  see  [13]. 

N'oca  chat  the  framework  described  here  can  a-iso 
handle  the  effects  of  multipath  propagation.  In  the 
presence  of  multipath  y^(z)  (and  N,(z))  will  have 
several  large  coefficients  corresponding  to  che  direct 
delay  and  che  multipath  dalaya.  Thus,  a  careful  examin¬ 
ation  of  the  coefficients  reveals  useful  information 
about  tha  multipath  structure.  Since  che  direct  path 
has  the  shortest  delay,  the  first  large  coefficient  of 
Ni(z)  will  correspond  to  the  direct  path  delay  and 
provide  T30A  information 

The  ARMA  modeling  approach  can  be  extended  to  a 
sultitarget  environment.  Here  che  system  consisting  of 
targets  and  sensors  will  be  represented  by  a  multi¬ 
output  (MIMO)  transfer  function.  A  simple  example  Is 
depicted  in  Figure  2.  The  equation  describing  the 
vector  of  measured  dace  y(k)  is  given  by 

y(z)  -  N(z)  X(x>  +  V(z>  (U) 

where  V(z)  represents  measurement  noise  and 

fz*TH  z-T2ll  , 

(or  a  ®ora  ganaral 

M(z)  -I  _T  _  -t__I  dalay  matrix) 


z-Tn  z’t: 


i?eCT!Wl  PROPAGATION 

RWl  WEL 


Figure  2.  The  Multitargec  Casa 
The  3ignal  vector  X(z)  Is  generated  by 

fd1(a)  0  ] 

X(z)  -  D_l(z)  U(z),  0(a)  -  j  .  (13) 

L  0  d2(l)J 

Tha  oucpuc  at  the  sensors  in  Che  multitargec  case 
can  be  written  as 

V(i)  -  N(z)  D(z)*1  U(z)  +  7(a)  .  (14) 

Noce  chac  as  in  the  single  target  case  the  numerator 
M(z)  contains  location  information  (TDOA's)  -while  the 
denominator  3(a)  contains  spectral  information.  If 
the  parameters  oi  N(z) ,  '(2)  are  given  the  aulcl- 
target  cracking  problem  is  basically  solved.  The 
question  is  how  to  estimate  these  parse  ters  from  the 
lata  •/( t)  in  a  unique  and  reliable  way.  As  nay  be 
expected,  the  question  of  uniqueness  is  somewhat  more 
complicated  than  in  the  single  target  case,  as  will  be 
discussed  next. 

3.  sncami  issi~s 

3.1  Tnlqueness 

The  identification  of  linear  systems  from  input/ 
oucpuc  data  requires  proper  parameterization  of  their 
transfer  function  (or  state-space  representation) .  In 
the  single-input  single-output  case  all  that  is  needed 
is  a  determination  of  Che  orders  of  the  numerator  and 
denominator  polynomials.  In  the  multivariable  case 
more  structual  information  is  needed  due  Co  the  fact 
that  the  numerator  and  .eliminator  polynomial  matrices 
can  be  transformed  in  various  nontrivial  ways  without 
changing  the  order  of  the  realization.  To  see  this, 
let  U(z)  denote  a  uniaodular  polynomial  matrix  [11]; 
then  H(z)«N(z)  Q(z)  -(S(z)U(z) )  (D(z)U(z)  T1 .  In 
other  words,  N(z)  *  N(z)CT(z}  and  5(z)  ■  0(z)U(z)  pro¬ 
vide  a  matrix  fraction  description  (MFD)  of  the  trans¬ 
fer  function  H(z),  of  the  same  order  as  N(z),D(z). 

To  obtain  a  unique  parameterization  of  the  transfer 
function  its  matrix  fraction  description  iN(z),  D(z>[- 
has  to  be  put  in  a  canonical  form.  Some  of  the  commonly 
•used  forma  are  the  Smith-McMillan  fora  and  tha  ?opov  or 
Polynomial-echelon  fora.  It  is  interesting  to  hots  that 
the  model  associated  with  the  muitltarget  cracking  pro¬ 
blem  (equations  (12)  (13))  has  a  diagonal  denominator 
matrix.  If  N(z),D(z)  are  irreducible,  this  implies 
uniquesness  of  the  MFD.  [To  see  this  consider  the  class 
of  uniaodular  matrices  U(z)  such  that  dlag{d^(z)}- 
u(z)  *  diagonal.  Only  J(z)  •  diag  iC^;  for  some  con¬ 
stants  will  satisfy  this  equation.  If  we  further 
require  chac  che  leading  coefficients  of  polynomials 
are  unity,  we  must  have  T(z)  -  I).  The  as.jsption  of 
irreduclbillty  of  N(z)  and  3(z)  will  be  generally 
satisfied,  since  it  is  unlikely  that  che  delay  structure 
N(z>  will  bear  any  relationship  to  the  target  signals 


D(z)  ■  Thus,  if  che  diagonal  structure  of  D(_z)  is  taken 
into  account  in  the  Identification  process,  and  if  the 
orders  of  Dtz) ,  N(z)  are  known,  a  unique  parameteriza¬ 
tion  of  che  transfer  function  will  be  achieved. 

For  the  tracking  problem  che  uniqueness  of  N(z), 
D(z)  is  not  too  important  as  long  as  the  transfer  func¬ 
tion  H (z)  can  be  uniquely  estimated.  It  is  possible, 
in  fact,  to  extract  the  information  needed  for  tracking 
directly  from  che  impulse  response  of  the  syscem  H(z)“ 
N(z)D(z)  "1 ,  -without  looking  at  the  numerator  and  denom¬ 
inator  polynomials.  3y  applying  an  impulse  to  che  input 
u.,  of  che  first  channel  (see  Figure  3)  we  will  obtain 
at  the  oucpucs  che  impulse  response  of  l/d-(z) 
delayed  by  various  amounts.  The  spectrum  of  target  No.  1 
can  be  evaluated  by  computing  the  power  spectrum  of  che 
impulse  response,  and  its  TDOA  by  cross-correlating  the 
two  outputs.  A  similar  procedure  can  be  carried  out 
for  Target  No.  2.  The  polnc  is  that  once  the  transfer 
function  from  u  to  y  is  known,  we  can  effectively 
separace  the  different  targecs! 


Figure  3.  Using  the  Model  to  Obcain 

Information  Abouc  Targec  No.  1. 

The  main  difficulty  in  the  tracking  problem  is 
that  we  have  to  estimate  H(z)  from  output  measurements 
only,  i.e.,  we  know  y(c)  but  not  u(t).  It  is  well- 
known  chat  for  stationary  gauasian  inputs  u(c)  the 
transfer  function  H(z)  can  be  determined  only  up  to  a 
spectral  equivalence  class.  In  other -words,  ail  -we  can 
say  is  that  H(z)  HU"1)*  -  S  (z),  the  spectrum  of  the 
data  y(t).  Making  che  usual  assumption  that  both  the 
poles  3nd  zeroes  of  H(z)  are  -within  the  unit  circle 
(i.e.,  a  stable,  minimum  phase  plant),  we  art  still 
left  with  nonuniqueness  of  che  spectral  factor  due  to 
the  existence  of  unit  spectral  aacrices.  Let  ^S£z) 
denote  a  polynomial  matrix  such  chac  S(z)  S(z  )  *  I. 

Then  H(z)  and  tl(z)S(z)  have  che  same  spectrum.  To 
set  that  these  exist  nontrival  unit-spectral  matrices 
let  Q  denote  an  orthogonal  matrix,  i.e.,  QQ*  • 
and  lac  R(z)  donee  a  diagonal  delay  matrix,  i.e., 

R(z)  •  diag  { x_fc3 )  for  some  integers  k, .  It  is  easy 
to  sae  that  arbitrary  products  of  matrices  of  chase  two 
types  are  a  unit  spectral  factor.  Consider  fo£  example 
S(z;-qR(z);  then  S(z)  ST(z)-QR(z)  R(z‘1)QT-QIQ‘-I. 

Tha  presence  of  unit  spectral  factors  of  the  ?.(z)  type 
has  the  intuitive  meaning  as  in  che  scalar  case;  the 
absolute  propagation  delays  can  noc  be  determined  from 
measurements  of  y(t)  alone.  Only  time  differences 
can  be  determined. 

In  order  to  be  able  to  separate  the  information 
related  to  the  different  targecs  it  is  essential  to  be 
able  to  get  rid  of  tha  nonuniqueness  introduced  b7  the 
unit  spectral  factor.  To  do  this  we  oust  use  the 
special  structure  of  che  fenooinator  matrix  3(z)  and 
assume  that  N(z),D(z)  are  irreducible  (i.e., 
minimal  realization).  First  note  that  N(z)D(z)  S(z) 
will,  in  general,  have  a  larger  order  than  N(z)T(z)-l 
unless  dec  3(z)*I.  Thus,  for  the  proper  choice  of 
oodel  order  we  must  have  S(z)«3.  In  ocher  words,  the 
only  nonuniqueness  left  if  we  restrict  our  attention  tc 


minimal  order  spectral  factors  of  S  (z),  is  generated 
by  orthogonal  craneformaslons .  Mexrnoce  that  H(z)  * 
S(*)0<z)'lq  •  S(z)D(z)(Q‘D(z)U(z))  ,  where  U(z)  is 

a  unimodular  matrix.  Since  we  'enow  chat  D(z)  Is 
diagonal,  che  nonuniqueness  Issue  reduces  to  the  ques¬ 
tion  of^vhac  class  of  matrices  L'(z)  and  Q  win 
obey  3‘  D(z)  U(z)  •  diagonal.  3y  the  uniqueness  of 
che  Smith-McMillan  form  [11]  It  follows  chat  che  only 
possible  modifications  of  D(z)  are  different  order¬ 
ings  of  the  polynomials  (d^(z) . d  (*)}•  In  ocher 

words  U(z)  -  Q  where  Q  Is  a  permutation  me  tlx  (It 
nay  also  include  arbitrary  sign  changes) .  This  fact 
has  che  following  Intuitive  interpretation:  since  the 
initial  labeling  of  che  targets  is  arbitrary,  we  can 
only  expect  to  define  H(z)  up  to  column  permutations. 
However,  once  a  given  labeling  was  determined,  H(z) 
will  be  fixed,  and  its  MFC  will  be  uniquely  determined. 

In  suamary:  using  the  diagonal  structure  of  D(z) 
and  assuming  Icnovledge  of  the  correct  (minimal)  order 
of  che  transfer  function  from  u(c)  to  y(c),  it  is 
possible  to  obtain  a  unique  (up  to  input  labeling) 
spectral  factor  of  Sv(z)  •  3(z)  H(z-1-).  The  impulse 
response  of  che  spectral  factor  H(z)  will  provide 
che  desired  cracking  parameters  (TDOA’s,  spectre). 

It  should  be  noted  chat  in  addition  to  the  structure 
of  D(z),  we  may  use  che  special  form  of  M(z)  which 
is  known  co  be  a  delay  matrix.  This  information  can  be 
used  co  relax  che  assumption  of  minimality  and  to 
obtain  uniqueness  for  more  general  problem  formula- 
cions  (see,  a.g. ,  Section  5) . 

3.2  Right  and  left  Matrix  Traction  Descriptions 


The  remsining  question  is  one  of  psraiMCer  estima¬ 
tion.  How  to  find  s  ssc  of  coefficients  N(z),  D(z) 
which  will  best  fit  (in  che  leaec-squares  or  maximum 
likelihood  sense)  a  given  cime-series  (y(e)K  Various 
techniques  have  been  proposed  in  che  literature  for 
estimating  che  parameters  of  ARMA  processes  of  che  form 
p  m 

y(t)  -  -  ^  At  y(t-i)  +  21  8t  •*<«-»  •  US) 

.Hose  of  che  work  in  this  area  seems  co  have  been  done 
in  che  context  of  system  Identification  (i.e.,  che 
known  Inpuc  ceae)  [13] -[201,  but  some  results  are 
available  on  time-series  modeling  (unknown  Input  case) 
[21].  'Jnfortunacely ,  moat  of  these  techniques  are  not 
applicable  co  the  estimation  of  H(z) ,  D(z).  Note  chat 
the  ARMA  model  In  equaclon  (13)  Is  naturally  related  co 
che  Laft  Matrix  Fraction  Description  (LMFD)  of  the 
transfer  function  from  u(t)  to  y(c),  while  H(z), 
D(z)  are  its  Right  Matrix  Fraction  Description 
(RMFD) .  To  see  this  we  rewrite  (13)  as 

A(z)  Y(z)  •  3(z)  U(z)  (16) 


A(z)  »  I  +  Ajz"1  ♦  ...  +  Ap*'P. 

3(z)  -  3„  *  S.z"1  +  ...  +  3  z*™ 
oi  a 

thus, 

H(z)  -  A(z)*1  B(z)  -  H(z)  DCz)*1  .  (17) 

One  way  of  estimating  H(z),  D(z)  Is  co  first  estimate 
A(z),  B(z)  using  any  of  the  techniques  described  in 
[15]-[21]  (see  also  Section  1),  and  then  compute  che 
RMFD  of  the  transfer  function  A(z)-)-B (z).  The  RMFD 
and  LMFD  are  releeed  by  A(z)M(z)  -  3(z)D(z)  ■  0. 
Writing  this  equation  in  terms  of  the  coefficients  of 
A.3.H  erd  D  leads  co  the  following  matrix  equation 
(written  for  simplicity  for  the  cast  m»p»2) : 


*1  Ao 

A2  *1  Ao 
A,  A1 

A, 


®2  B1  3 


where  •  I.  Mote  that  if  we  fix  the  leading  coef¬ 
ficient  of  D(z) ,  say  D  "I,  equation  (18)  can  be 
rewritten  as  a  matrix  equaclon  with  number  of  unknowns 
equal  co  the  number  of  equations,  and  can  therefore  be 
solved  for  {M  ,D.}.  This  matrix  equation  involved  a 
(modified)Sylvesier  matrix  with  a  very  special  structure. 
Using  this  structure  Rung  at  el.  [ 22 ] — [ 24 ]  developed 
efficient  algorithms  for  going  from  LMFD  co  RMFD  and 
vice-verse.  These  algorithms  ucilize  orthogonal  trans¬ 
formations  and  are  claimed  co  have  good  numerical  pro¬ 
perties.  The  problem  of  changing  right  to  left  matrix 
fraction  descripcione  is  closely  releeed  co  finding  che 
greetesc  coenon  divisor  of  polynomial  matrices  end  co 
computing  e  minimal  basis  for  che  space  spanned  by 
such  macrlces  [11]. 

A  preferable  solution  would  be  to  estimate  directly 
the  parameters  (N  ,D  }  of  the  RMFD.  This  can  be  done, 
for  example,  by  using  a  maximum  likelihood  or  prediction 
error  formulation  which  leads  Co  e  nonlinear  optimiza¬ 
tion  problem.  Such  a  aoluclon  will  in  general  be  com¬ 
putationally  expensive.  Some  preliminary  Investigation 
seems  co  indicate  che  possibility  of  obtaining  estima¬ 
tion  algorithms  for  N(z),D(z)  which  are  of  a  compar¬ 
able  level  of  complexity  Co  algorithms  for  estimating 
A(z) ,3(z) .  However,  chis  investigation  is  not  yet 
complete. 

It  should  be  recalled  chat  the  desired  target 
information  eta  be  obtained  directly  from  the  impulse 
response  of  che  system,  which  can  be  computed  either 
from  the  LMFD  (H(z)-A(z)-l-3(z))  or  from  che  RMFD 
H(z)»N(z)D(z)"1  as  Indicated  in  Section  3.1.  The  mein 
reason  for  wanting  Co  esclmate  che  RMFD  is  Chat  ic 
displays  claarly  che  special  structure  of  che  system 
which  Is  nesded  co  fores  a  unique  solution  of  the  spec¬ 
tral  factorization  problem.  If  che  special  structure 
of  H(z) ,D(z)  could  be  mapped  into  an  easily  specified 
special  structure  for  A(z),3(z),  we  could  work  directly 
in  terms  of  Che  LMFD.  Unfortunately ,  che  structure  of 
A(z),B(z)  Induced  by  D(z)  being  diagonal  and  M(z) 
being  e  delay  matrix,  seams  to  ba  complicated  and  dif¬ 
ficult  to  use.  Ic  should  also  be  noted  chat  it  is 
desirable  co  enforce  the  special  structure  already  in 
the  estimation  process.  This  can  be  done  easily  in  the 
RMFD,  buc  not  in  the  LMFD,  which  makes  It  even  more 
important  co  develop  s  direct  estimation  technique  for 
N(z),D(s). 

In  suaoary,  several  approaches  to  che  estimation 
of  cht  modal  paramsters  wars  discussed:  (i)  Direct 
estimation  of  chs  RMFD,  with  D(z)  forcsd  co  bs  diag¬ 
onal,  (li)  Estimation  of  Chs  LMFD  using  availabla  ARMA 
modeling  techniques.  Transformation  of  che  esc.^eted 
LMFD  into  RMFD  using  available  techniques.  Checking 
the  structure  of  D(z)  and  finding  the  transformation 
nssdad  to  maka  it  diagonal.  This  transformation  will 
pin  down  s  unique  epeetrel  factor,  (ill)  Estimate 
A(z),3(z)  with  their  form  forced  co  have  the  particular 
structure  Induced  by  K(z),D(z).  Then  compute  che 
impulse  response  of  H(z)  •  A_i(z)S(z)  and  from  it  the 
TDOA’s  end  cargee  spectra.  Of  these  approaches  only 
che  second  is  fully  developed  at  che  current  time,  ee 
described  in  Section  4. 


3,3  Innovation;  Rtortsantatloa 

In  eh*  dt*cu*»loa  to  Ur  w*  hav*  not  mentioned 
explicitly  eh*  nuabar  of  eargecs  or  t«ntort.  Lae  ua 
aaauaa  ehac  ch«r*  ar*  N_  targets  tad  .v.  tanaors 
and  eh*e  3,  N*.  Consider  eh*  tpaecru a  of  y(e) 

given  by  aquatint  (Li) : 

S  (a)  -  X(z)0(z>*1  afr'L)*T  .V(z*ST  +  Z  (19) 

y  v 

whara  Z  danoeaa  eh*  covarlanea  matrix  of  eha  measure¬ 
ment  not**  v(c).  Note  chae  if  Iv*0  aad  st<Ns, 

3  (a)  will  aoc  b*  a  full  rank  matrix!  (An  equivalent 
way  of  scaclng  chi*  face  is  eo  tay  chae  eh*  inaovaeioo* 
representation  of  y(e)  will  b*  singular,  aad  son*  of 
eh*  components  of  eh*  ioaovacioaa  vaeeor  e(e)  will 
b*  linaarly  dependent) .  In  fact,  eh*  rank  of  eh*  covar- 
ianct  matrix  5  (z)  (or  eh*  innovations  covarianc* 

aaerix  R")  will  b*  equal  to  en*  nuabar  of  carg *ea,  aa 
can  b*  saan  froo 

Sy(t>  •  S(!)0(:)’1  SCi*1)1  (10) 


Thu*,  eh*  nuabar  of  eargecs  can  b*  aaclaatad  by  eaaeiag 
Ch*  rank  of  5y(z)  (R1) .  Thia  caat  will  ba,  of  court*, 
sensitive  tine*  any  measurement  nolaa  will  aaka  Z  >0 
and  caua*  S  (z)  (R~)  eo  bacoaa  full  rank.  However,  by 
looking  ac  th*  ralativ*  magnitude*  of  eh*  tig«nvalu*t 
of  eha  covarianc*  macrix,  w*  can  obtain  a  aora  robuac 
estimate  of  eh*  nuabar  of  cargoes.  This  procedure  la 
daacrlbad  in  aora  detail  in  Saccion  1.  In  g«a*ral,  wa 
will  work  with  a  "squara"  eransfar  funecioa,  i.*., 

N(z)  and  3(z)  (or  A(z)  and  3(z))  art  MjXMj  aatrieat, 
and  dacid*  on  th*  actual  nuabar  of  cargoes  afear  per¬ 
forming  th*  nodal  fitting.  If  w*  nay  think 

of  this  as  having  actual  targets  and  Jfj-W, 

ficcicious  eargaca  which  ganaraea  no  signals  (its., 
th*  corrssponding  inputs  u^O-O) . 


F(z)«N(z)  ^YC*)  and  spectral  fiicaring  U(z)»D(z)F(z)  , 
performed  sinuleanaoualy  for  all  targets! 


A.  THE  ALGORITHM 

A  candidate  algorithm  for  estimating  th*  para¬ 
meters  A^.3^  2*  described  below.  This  algorithm 

it  the  multivariable  version  of  eh*  Recursive  Maximum 
Likelihood  (RML)  cechlqua,  iatcribad  in  [25]  for  ch* 
tingle  inpuc/oucpue  cat*. 


Lac  us  rewrite  eh*  ASHA  modal  of  aquadon  (15)  aa 


(21) 


where 


i  -  [A . A  B . BJ, 

1  p  O  S 

an  NjX(p4myl)N^  aaerix 


t  -V  !• 


-y. 


'e  *  'c-1’  'c-p’  ~t-l . t-»-l 

an  1  x  (p-ran-l)Sj  vector 


or  it  can  b*  decoupled  into  Ng  equations  of  eh*  eyp* 

't\  ’  5c  +  s't  ’  3-1 . (22) 

where  d3  is  eh*  J-eh  column  of  3  . 

The  estimation  algorithm  for  i  is  given  by 

*  -i  *  xtti  £c+i*  J*1 . <a) 

where  is  a  gain  vector  (common  to  all  of  ch* 

d3  ’s) ,  and 

i+1  •  rc+l  -  Je+l  3-1 . Ns  <24> 


Th*s*  aquation*  can  b*  written  in  a  more  compact 
form  as 


3.1  3«amforaing  Incarpraeaeion 

A  standard  approach  to  passive  tracking  ia  baam- 
formlng.  In  ica  simplest  form,  beamforming  consists 
of  forming  a  sum  of  delayed  sensor  outputs,  to  gat  a 
scalar  signal  f(c)  which  chan  undargoa*  furthar  pro¬ 
cessing.  In  ch*  single  target  caae,  this  can  ba  viewed 
as  forming  an  inner  product  of  cha  data  vector  Y(z) 
wich  a  "steering"  vector  W(z), 

F(z)  •  W(z)T  Y(z)  •  «(z)T  (31(z)  XU)  *  V(z)> 

Th*  steering  vector  it  chosen  to  that  it  "remove*''  eh* 
numerator  polynomial  which  rapraaants  eh*  propagation 
modal,  1.*.,  chooa*  V(z)  so  chat 

W(z)T  N(z)  •  z*T 

and  chan  th*  received  signal  5(c)  la  juse  a  dalaytd 
varslon  of  eh*  targat  signal  x(c).  To  saa  thia  more 
clearly,  conaldar  »(*)  aa  given  in  aquation  (6),  and 
lac 

•d(z)T  -  (nz'T:  z*Tl]/2n  . 

Than  clearly  W(z)T  N(z)  ■  z'^Tl+r23.  ia  th*  multl- 
eargat  caa*  X(z)  will  b*  a  aaerix  which  is  chosen  so 
that  V(z)T  S(t)  •  diag  (z  Tl}.  Thua,  eh*  operation  of 
baamformlng  can  b*  thought  of  aa  finding  th*  invars*  of 
N(z)  in  a  generalized  Sanaa.  Moca  chat  baamformlng  is 
vary  cloaaly  raaltad  to  ch*  RMFD  representation.  In 
fact,  th#  whitening  filter  UU)»C(z)SU*l)Y(*)  can  b* 
considered  as  a  combination  of  baamformlng 


Vl  ■  *c  *  *fl  £t+l 

(25a) 

Ec*l  ■  yc+l  Vl 

(25b) 

Th*  gain  vector  K  can  b*  computed  by 

WptW(V*L.l?cW  ’  ?c+i’c+i 

(26a) 

vr^t-'tW^t^  vlipt'vi,1At 

(26b) 

whera  \t  is  an  exponential  weighting  factor. 

and 

i 

t 

la  a  filtered  version  of  th*  *  vector,  1.*., 

y.e 

la 

th*  p  vectors  are  replaced  by  pr«-fllt*r*d  v«r*ion* 

y  and  s  where 

7t  *  D_1(z)  yt  ,  et  -  D*1(s)  £t, 

(27a) 

DU)  •  1  +  01z”1  ♦  *  0kz_ic 

(27b) 

Tha  pre-filter  D(z)  i*  typically  chosen  a*  D(z)« 
B(x).  Howavar,  wa  found  that  Improved  convergence  can 
somatimas  b*  achieved  by  using  0(z)«A (z/a) ,  i<  1. 

Saa  [13]  for  more  detail*.  ~ 

Th*  innovations  sequence  produced  by  cha  RML 
will  generally  hav*  eorralacad  ancrlas  (i.*.,  a  noo- 
dlagonal  innovations  aaerix).  Sine*  w*  aaaum*  ch*  car- 
gat  signals  to  ba  uncorralatad,  it  ii  nacaaaary  to  work 
with  ch*  normlalzad  innovations  sequence  -- 


(28) 


where  S  ”  is  tha  (lower  triangular)  squats  root  of 
tha  innovation  covariance.  R'  can  be  computed  recur¬ 
sively  by 


*  a_Mc)  *c  +  wt  sc  Jt 

“c«.  -  We/(WC  *  V  (29b) 

or  in  square  root  form 

f-^C^2]  orthogonal  _ 

—  ,T  J  transformation  L  ^  j 

tha  normalized  innovations  will  ba  coapucad  by  solving 

K,Z  K  -  S  ‘31> 

which  can  ba  dona  by  back  substitution  sines  R  ^ “  is 
upper  triangular.  Spaclal  cars  has  to  ba  taken,  since 
R'  say  ba  singular,  i.e. ,  sous  of  tha  diagonal  ele¬ 
ments  of  R'/_  may  be  tero  (or  very  small).  In  that 
case,  Che  corresponding  elements  of  : £  will  ba  sac 
co  taro,  this  can  ba  dona  easily  in  che  back  substi¬ 
tution  routine,  tha  auaber  of  diagonal  eleaeocs  of 
R_/‘  which  are  graacar  chan  some  threshold  value, 
will  lndlcace  tha  number  of  cargacs,  as  was  explained 
earlier. 


to  eacimata  TDOA's  and  target  spectra,  tha  impulse 
response  of  H(z)«A(z)~  li(z)  R"  will  ba  coapucad 
for  each  inpuc/oucpuc  pair.  Sy  cross-correlating  the 
impulse  responses  obtained  for  a  given  input  i,  the 
TDOA's  corresponding  co  target  auaber  1  will  ba 
obcained.  tha  power  spectrum  of  the  impulse  response 
will  provide  an  estimate  of  tha  target  spectrum,  there 
ia,  however,  one  remaining  problem:  the  transfer  func¬ 
tion  H(z)  may  ''scramble  up"  the  targets  by  an 
orthogonal  transformation.  In  other  words,  the  correct 
transfer  function  will  generally  be  H(z)Q,  where  Q  is 
an  orthogonal  transformation  that  needs  to  be  deter¬ 
mined  (alternatively:  fin^,  che  correct  way  of  caking 
che  matrix  square-root  R~  ")• 


In  the  case  where  is  small  it  is  possible  to 

taka  a  general  orthogonal  transformation  matrix  and 
look  at  the  cross-correlations  of  che  Impulse  reponses, 
obcained  for  different  rotations.  If  che  wrong  rota¬ 
tion  is  used,  che  cross -correlation  function  will  have 
multiple  peake  correepondlng  to  several  cargacs.  If 
che  correct  orthogonal  transformation  is  used  only  a 
slngls  peak  will  occur  corresponding  co  a  slngla  tar¬ 
get.  this  is  somewhat  similar  co  beam  steering 
(except  that  it  Is  done  in  e  different  space) . 


A  xare  systematic  way  for  decetalnlag  tha  ortho¬ 
gonal  transformation  Q  is  to  compuce  the  RMFD  of 
a(z)'4(z)R'/‘.  One  technique  for  doing  this  la  pre¬ 
sented  in  (221  and  will  not  be  repeated  here.  Once  the 
RMFD  h(z),9(z)  is  calculated,  Q  can  ba  determined 
as  che  craaaofzmeclon  needed  co  make  D(z)  diagonal 
(or  as  nearly  diagonal  as  possible) . 


Finally,  we  should  nots  that  che  algorithm  pre¬ 
sented  here  was  chosen  mainly  for  convenience  of  Imple¬ 
mentation.  Many  ocher  algorithms  can  be  used  to  sati¬ 
nets  and  it  Is  not  deer  at  this  time  which 

one  is  best  for  multitargat  tracking  applications . 
However,  che  RKL  ttchniqua  saama  adequata  for  casting 
the  proposed  modeling  approach.  Some  preliminary 
simulation  results  are  given  in  the  Appendix. 


5.  SOME  EXTENSIONS 

the  multitargat  cracking  problan  was  formulated 
in  Sactlon  2  for  flxsd  cargacs  with  AS  spectra,  and  a 
simple  propagation  modal.  In  this  section  wa  briefly 
describe  how  che  ASMA  modeling  approach  can  be  ascended 
to  ocher  daaaes  of  cracking  problems. 


When  cargacs  and  sensor  locations  are  fixed  in 
space,  tha  modal  presented  in  Section  2  will  ba  time- 
lavarlanc.  (This  follows  from  eba  assumptions  of  the 
stationary  target  signals  and  a  time- invariant  propaga¬ 
tion  mediia) .  Raise iva  target  motion  will  have  two 
important  affects:  (1)  Tha  system  relating  y(t)  to 
x(e)  will  be  time  varying  and  (11)  tha  signals  gen- 
eraced  by  the  cargee  will  experience  doppler  shifts, 
tha  first  affect  means  that  we  now  heve  co  identify  an 
A8MA  modal  with  time-varying  coafflciaats.  (In  face, 
only  the  N  coefficients  are  elms  varying) .  This 

is  a  difficult  task  even  in  cha  single  channel  case, 
although  it  has  bean  dona  successfully  in  various 
applications  such  as  speech  processing.  In  many  situa¬ 
tions  (e.g.,  sonar  systems)  target  motion  ia  sufficiently 
slow  compered  co  che  sampling  rate,  to  thee  che  para¬ 
meter  estimation  algorithm  will  be  able  co  creek  para¬ 
meter  changes.  More  precisely,  ee  long  as  cha  target/ 
sensor  geometry  can  ba  considered  fixed  over  che 
time  interval  needed  to  get  a  reasonable  location  esti¬ 
mate  (integration  time)  techniques  which  work  wall  for 
fixed  cargacs  will  3 till  be  expected  to  work.  If  the 
geometry  changes  sufficiently  fast,  a  time-varying  pro¬ 
blem  formulation  is  unavoidable.  For  example,  we  may 
cry  co  parameterize  the  coefficient^  S  of  che  delay 
matrix  by  S. (t)  -  N,  +M,  ,t+M,  ,t“.  The  coefficients 

a  If*  If* 

N.  M,  ,  S,  ,  can  be  related  to  the  target  location 
10,  1,1,  i, 2 
and  its  valocity. 


In  many  applications,  the  signals  rscelved  from 
moving  targets  can  undergo  significant  doppler  shifts, 
even  if  cha  carget/sensor  geometry  is  only  slowly  time- 
varying.  These  frequency  shifts  provide  information 
about  target  velocity  and  about  its  location  (e.g.,  by 
looking  at  differences  of  doppler  shifes  between  pairs 
of  sensors),  thus,  1:  is  important  co  see  how  doppler 
effects  fit  into  the  modeling  framework.  Consider  an 
AR  process  x(t)  generated  by  a  target,  where 
x(t)"(l/d(z))  u(c).  For  narrowband  processes,  it  can 
bs  shown  thee  the  doppler  shifted  signal  will  still  be 
an  AR  process,  except  that  the  AR  coefficients  d^ 
era  changed.  In  feet.  If  d(z)«(z-z1) . . . (z-z  )  then 
its  doppler  shifeed  version  will  be  d(z)»(z-z^; . . .  (z-zj) 
where  Cha  parameter  a  Is  determined  by  cergec 
velocity  tod  propagation  velocity.  Thus ,  the  model 
for  the  observed  signal  will  heve  che  form 


y(i)  -  H(z)  0(z)  +  V(z)  (32s) 

H(z)  -  [z_TlJ/dtJ(z)J,  Ks  x  ST  matrix  (32b) 

where  d  (a)  is  ths  AR  representation  of  che  signal 
generac«dJby  target  1,  ea  aaan  by  sansor  J,  and  T^ 
is  tha  daisy  from  target  1  to  sensor  J.  Note  that 
che  parameter  of  d^. (z),  )*1,  ....  M.  are  not 
Independent,  but  are  related  in  e  special  way.  An 
examination  of  these  parameters  provides  not  only  a 
spectral  estimate,  but  also  information  about  cergec 
velocity  and  location,  the  structure  of  this  modal 
has  many  interesting  features  which  are  still  under 
investigation.  Of  particular  interest  is  cha  question 
of  how  to  obtain  a  unique  spectral  factor,  la  view  of 
cha  feet  that  D(z)  is  no  longer  diagonal.  As  may  be 
expected,  uniqueness  can  be  achieved  under  fairly  mild 
conditions.  Estimating  cha  parameters  of  such  modais 
from  noisy  data,  in  a  way  that  takas  into  account  chair 
spaclal  structure  is  not  vet  fully  developed. 


I 


5.2  ABft  Source* 

While  AR  models  provide  aa  adequate  signal  repre- 
saacacion  la  a  larga  elaaa  of  practical  eases,  lc  la 
sons  class  necessary  to  consider  the  nor a  general  ARMA 
representation.  This  is  particularly  true  for  wide¬ 
band  processes  or  narrowband  processes  la  additive 
noise.  The  signal  aodel  will  have  In  chls  case  the 
fora 

X(t)  •  5(z)  Sft)"1  U(z)  03) 

where  2t( c) ,  5(z)  are  both  diagonal, 

and  . 

Y(z)  -  (»(*)  K(z) )  Sfr)'1  tJ(a)  +  V(z).  (3A) 

Note  thac  the  auaarator  of  the  RMTD  la  now  a  aixture 
of  che  delay  aacrix  and  Che  spectral  zeroes  of  Che 
cargecs.  however,  SCz)  Is  still  diagonal  and  chare- 
fore  nost  of  che  discussion  la  Section  3  still  holds 
true.  The  TDOA's  and  cargec  spectra  can  still  be 
esclaaced  from  che  lapulse  response  of  che  estiaaead 
aodel.  3ecauae  of  che  diagonal  structure  of  N(z) 

It  Is  theoretically  possible  co  separata  S(z)  and 
si(z)  _by  looking  for  common  factors  la  the  coluans  of 
3(z)  ,V(z>.  In  sut ranary ,  having  an  ARMA  spectral 
aodel  does  not  complicate  significantly  the  aodeliag 
approach. 


ARMA  Propagation  Models 


The  siaplesc  fora  of  aulclpach  propagation  is 
caused  by  a  signal  being  reflected  by  tone  object 
which  lies  outside  the  direct  line  of  propagation, 
causing  aa  exersneous  delayed  version  of  che  signal 
co  arrive  at  che  sensors.  This  type  of  aultlpeth  pro¬ 
pagation  can  be  represented  by  a  MA  aodel  of  Che  type 
used  in  earlier  sections.  Often,  a  aora  co^>llcated 


type  of  propagation  occurs.  The  signal  can  undergo 
aulciple  reflection*  of  che  type  encountered,  for 
instance,  whan  a  sound  wave  propagates  In  a  room  (a.g. , 
che  "cocktail  party"  problem:  locating  people  who  are 
speaking  inside  a  room).  The  sound  is  raflscced  from 


one  wall,  bounces  off  che  opposite  well  and  again 
from  the  first  wall.  Similar  sffacts  occur  when 


seismic  waves  propsgats  in  the  sarth  or  when  sound 
propagates  in  the  ocean  (reflected  from  the  waeer/air 
and  waccr/oesan  floor  interfaces) .  This  phenomena  is 


best  modeled  by  an  AR  cype  aodel. 


In  general,  a  reallatic  propagation  aodel  will 
combine  both  types  of  multipath  propegeclon  and  will 
therefore  be  an  ARMA  aodel,  i.e., 

T(t)  -  N(z)  D(z)'1  X(z)  +  V(z)  - 


.  S(z)  [D(z)  D(z)]'1  U(z)  +•  V(z)  (35) 


where  X(z) ,  D(z)  are  the  propagation  aodel  and 
5(z)  is  che  signal  aodel. 


The  structure  of  the  propagation  aodel  depends, 
of  course,  on  che  particular  application  being  con¬ 
sidered  and  on  the  physical  properties  of  the  propaga¬ 
tion  medium.  In  its  simplest  fora  D(z)  will  have  a 
diagonal  fora,  aaaning  that  che  propagation  from  a 
given  cargec  co  all  sensors  has  the  seme  kind  of  pro¬ 
pagation  effaces,  except  for  different  delays  which 
are  represented  in  M(z).  In  general,  aora  complicated 
forae  of  propagation  are  possible  in  which  case  3(z) 
will  be  non-diagonal. 


When  3(2)  is  diagonal,  che  structure  of  the 
transfer  function  H(r)  •  S(z)  (5 (t)  D(z)))  1  is 
similar  to  the  pure  delay  case.  The  only  problem  is 
that  the  signal  spectrum  is  mixed  up  with  the  poles 


of  che  propegeclon  medium.  Under  certain  situations 
lc  aay  ba  possible  to  sepsrat#  the  two.  Consider  the 
situation  where  the  time  dependence  of  the  signal 
spectrum  and  the  propagation  model  are  diffsrsnt.  For 
example,  consider  a  none cations ry  source,  like  speech, 
which  stays  fixed  in  specs  so  chat  cho  propagation 
aodel  is  cima-lnvarianc.  By  repeating  che  computation 
of  5(z)  3(z)  over  several  time  intervals  and  comput¬ 
ing  its  roots,  3(z)  will  be  found  from  chose  roots 
which  do  not  change,  while  3(z)  will  correspond  to 
che  roots  which  era  varying  from  one  time  interval  to 
another.  As  another  example,  consider  a  moving  cerget 
thee  causes  the  propagation  aodel  to  change.  Here 
5(i)  will  have  fixed  roocs  and  D£t)  roocs  which  vary 
over  time.  If  either  D(z)  or  D(z)  ere  nondlagonal 
the  structure  of  the  cransfer  problem  becomes  aora 
coepliceted  and  has  to  be  examined  more  carefully. 

This  will  not  ba  done  hare. 

5.4  Olrect  Estimation  of  Source  Location 

The  tee  of  TDOA’s  contain  all  the  information 
required  co  locate  the  targets,  i.e.,  to  find  their 
bearing/ range,  location  estimation  is  typically  par- 
formed  in  cwo  steps:  First,  estimate  che  TDOa’s  (2  ^ ) 
using  the  algorithm  described  earlier,  or  e  more 
standard  generalized  cross-correlation  technique.  Then 
compute  target  locations  from  the  sat  of  , ’s  using 
che  approach  presented  in  [12],  or  the  hyperbolic  llne- 
of-posltion  technique. 

A  different  approach  will  be  to  cry  and  estimate 
directly  cha  target  coordinates.  Since  che  functional 
relationship  between  the  target  coordinates  and  the 
TDOA’s  Is  known,  it  is  possible  co  re-foraulaco  che 
problem  as  a  coordinate  estimation  problem.  Some  pre¬ 
liminary  studies  seem  co  indicate  the  feasibility  of 
the  direct  approach,  but  no  conlcuslve  resulcs  ere 
available  at  this  time. 

6.  CONCLUSIONS 

We  introduced  e  framework  which  relates  various 
problems  thac  arise  in  multltarget  tracking  Co  Che 
perazMtera  and  the  structure  of  e  related  multivariable 
linear  system.  This  approach  activates  some  interesting 
questions  regerdlng  the  modeling  of  vector  ARMA  pro¬ 
cesses  under  various  structural  constraints.  The 
results  presented  here  ere  only  preliminary  end  che 
objective  was  mainly  co  lay  the  groundwork  for  further 
research . 

The  proposed  approach  has  eavarel  features  which 
are  Important  from  a  practical  standpoint.  Including: 
creating  multiple  targets  without  performing  a  eaparace 
association  problem  (targtt  assocatlon  as  well  as  line 
association) ,  simultaneous  estimation  of  TDOA  end 
target  spectre  and  possibility  of  handling  multipath 
effects.  However,  considerable  tea ting  and  performance 
evaluation  needs  co  be  done  before  the  usefulness  of 
the  ARMA  modeling  approach  can  be  fully  realized. 
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APPENDIX  C 


MTS  PARAMETER  ESTIMATION  ALGORITHMS 


A  number  of  algorithms  for  estimating  the  parameters  of  AR,  ARX  and  ARMA 
models  were  coded  and  tested  during  the  MTS  project.  These  algorithms  can  be 
divided  into  three  classes:  recursive  parameter  estimation  algorithms  of  the 
RML  type,  non-recursive  AR  and  ARMA  algorithms  and  recursive  lattice 
algorithms.  In  this  appendix  we  list  the  algorithms  that  are  currently 
available  and  either  give  a  brief  description  or  a  reference  to  a  more 
complete  description. 

C.l  RECURSIVE  IDENTIFICATION  ALGORITHMS  (IDNI) 

1.  Recursive  Least-Squares  (RLS) 

This  algorithm  estimates  the  parameters  of  the  ARX  model 


h  '  £ai  Vi  *£bf  Vi  *  \ 


where  ut  is  a  white  noise  process.  The  algorithm  is  given  by 
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2.  Extended  Least  Squares  (RM1) 

This  is  the  RM  algorithm  described  by  equations  (2)  -  (7)  with 
D^(  z)  3  1  and  with  e^  replaced  by  3  yt  -  <>^9^. 

3.  Extended  Least  Squares  (RMLP) 

This  is  the  RM  algorithm  described  by  equations  (2)  -  (7)  with 
DtU)  ■  1. 

4.  Recursive  Maximum  Likelihood  (RM) 

This  is  the  RML  algorithm  described  by  equations  (2)  -  (7). 

5.  Square-Root  form  of  the  RM  (NORM) 

This  is  the  algorithm  described  in  section  2.1  (11  i). 

6.  Stochastic  Approximation  RLS  (LMS) 

This  is  a  version  of  the  Widrow-Hoff  gradient  research  LMS  algorithm. 
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C.2  NONRECURSIVE  ALGORITHMS  (STAN) 

1.  Maximum  Entropy  Method  (MEM) 

This  Is  an  implementation  of  Burg's  technique  for  autoregressive  spectral 
estimation.  The  code  used  was  taken  from  [21]. 
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2.  Covariance  Method  (COV) 

This  is  an  implementation  of  the  covariance  method  of  speech  analysis  as 
given  by  [22,  pp.  5£].  The  resulting  AR  coefficients  are  used  to  compute 
the  spectrum. 

3.  Modified  Yule-Walker  Method  (MY VO 

This  is  an  algorithm  for  estimating  ARMA  parameters.  The  algorithms 
estimates  first  the  AR  parameters,  using  the  following  modified  Yule- 
Walker  equations: 
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This  overdetermined  system  of  equations  is  solved  by  a  least-squares  algorithm 
involoving  singular  value  decomposition  (ILSQF  in  the  IMSL). 

The  MA  coefficients  are  confuted  next  using  the  following  equation 
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The  spectrum  is  finally  confuted  by 


cM  .  b-NAzNA  +  +  b0  +  b-lz  1  bNAz'NA 

- 7T - ttt - — - mr 

(1  +  aj  z  ♦...+  aNAz  )  (1+8^2+...+  aNAzN*) 

where  z  •  exp 
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C3.  RECURSIVE  LATTICE  ALGORITHMS  (LADD) 


1.  The  Normalized  AR  Lattice  (ARN) 

This  algorithm  Implements  the  normalized  algorithm  described  in  [23].  The 
AR  coefficients  are  computed  by  Levinson's  algorithm  rather  than  by  the 
true  whitening  lattice.  See  [24]  for  details. 

2.  The  Unnormalized  AR  Lattice  (ARU) 

This  is  the  normalized  lattice  with  its  output  properly  scaled  to  give  the 
unnormalized  prediction  error.  See  [24]  for  details. 

3.  The  Two  Channel  Lattice  {KNWNl 

This  is  the  ARMA  lattice  for  the  known  input  case  [25]. 

4.  The  ARMA  Lattice  (ARMA) 

This  is  the  ARMA  lattice  for  the  unknown  input  case  with  the  prediction 
error  fed  back  to  the  input  [25]. 

5.  The  AM.  Lattice  (AM.) 

This  is  the  ARMA  lattice  for  the  unknown  input  case  with  the  residual  fed 
back  to  the  input  [25]. 

6.  The  Normalized  Lattice  (ARL) 

This  is  the  normalized  AR  lattice  with  the  true  whitening  filter  for 
computing  AR  coefficients  [26]. 

7.  The  Joint  Process  Lattice  ( ARJ ) 

This  is  the  normalized  joint  process  AR  lattice.  The  AR  coefficients  are 
computed  by  a  true  whitening  filter  [26]. 

8.  The  Sliding  Window  Lattice  (ARS) 

This  is  the  sliding  window  AR  lattice  form  [26]. 
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APPENDIX  D 


TRI ANGULAR I ZAT ION  BY  HOUSEHOLDER  TRANSFORMATION 


The  following  routine  will  take  an  LxN  matrix  and  upper- triangul arize  its 
first  M  columns. 


Note:  M  <  min  {N,L> 

Input:  an  LxN  matrix  S(I.J) 

M  =  #  of  columns  to  be  triangul ari zed 

Output:  an  LxN  (partially)  upper  triangular  matrix  S(I,J) 

DO  4  J*l,  M 

o  *  Z  @  Z  =  zero 

DO  1  I*J,  L 

v ( I )  *  S(I,J)  9  v(N)  is  a  work  vector 

S(I,J)  *  Z 

1  o*o  +  v{  I )  **  2 
IF  (a  <  Z)  GO  TO  4 

Comment:  If  o*Z,  column  J  is  zero  and  this  step  of  the  reduction  is 
omitted. 

o  *  SORT  la) 

IF  < v(J)<Z)  o  *  -o 
S(J,J)  *  o 


-  46  - 


0  ONE  *  1. 


v(J)  3  v(J)  -  a 
a  3  ONE/ (a  *  v(J)) 

00  3  K=J+1,  N 
a  3  Z 

00  2  W,  L 

2  a  3  a  +  S(I,K)  *  v<I) 
a  *  a  *  a 

DO  3  I3J ,  L 

3  S( I,K)  3  S(I,K)  +  a  *  v(I) 

4  CONTINUE 

In  the  system  identification  problem,  we  want  full  triangularization  of  a  square 
matrix,  thus  L=M=N. 
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